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ABSTRACT 


Intraframe  image  data  compression  systems  are  analyzed  In  tills 
thesis  using  stochastic  modeling  concepts.  The  formulations  of 
stochastic  image  models  are  obtained  from  different  classes  of 
partial  differential  equations.  Their  application  to  the  coding 
problem  shows  the  c  unectlon  between  predictive,  hybrid  and  transform 
coding  schemes.  The  resulting  coding  schemes  are  evaluated  in  terms 
of  tangible  system  terms  such  as  signal  to  noise  ratio,  mean  square 
error  and  rate  distortion  curve. 

Most  of  the  existed  intraframe  coding  schemes  are  either  ad  hoc 
or  lack  a  complete  mathematical  basis.  By  paying  attention  to  the 
boundary  points  of  an  Image  random  field,  it  is  shown  that  the 
remaining  image  random  field  can  be  fairly  described  by  stochastic 
image  models  which  in  turn  can  be  processed  via  well-defined  fast 
transform  algorithms  (Sine  and  Cosine  transforms  are  used  here). 

The  eigenvalues  associated  with  these  transform  algorithms  are 
important  for  coding  application.  To  properly  handle  a  nonstationary 
Image  random  field,  a  new  method  based  on  representing  an  image  as  a 
composition  of  two  sources,  viz.,  stochastic  and  deterministic,  is 
suggested.  The  stochastic  part  Is  dealt  with  by  known  statistical 
properties  of  stochastic  Image  monels  and  the  deterministic  part 
representing  features  such  as  edges,  is  obtained  as  the  residual 
after  the  stochastic  part  Is  removed.  Two  simple  adaptive  schemes, 
t.e..  Adaptive  Variance  Estimation  and  Adaptive  Classification,  are 
considered  for  updating  the  model  parsm-ters  with  the  variations  in 
image  statistics.  The  computational  complexity  is  Increased  only 
marginally,  but  the  Improvement  in  the  reconstructed  image  quality 

vii 


is  substantial. 


J 


Simulation  studies  are  performed  on  three  image  data  conslntlng 
of  Ctrl,  Earth  and  Chemical  Plant  Images.  The  Results  indicate  that 
the  proposed  schemes  can  be  practically  Implemented  and  will  reproduce 
a  good  iaMge  at  cosiparable  signal  to  nolae  ratio.  Hybrid  coding  of 
noisy  image  and  sensitivity  of  human  viewer  to  the  sharp  features 
are  also  partially  addressed. 
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CHAPTER  ONE 
INTRODUCTION 

1*1  The  Problem 

Digital  inage  processing  has  attracted  attention  and  found  applica¬ 
tions  In  diverse  fields  such  as  transmission  of  satellite  pictures, 
enhancement  of  biomedical  images,  radar,  seismic  signal  communication  and 
many  others.  Cenerally  speaking,  digital  image  processing  deals  with 
two-dimensional  data  which  is  first  sampled  in  spatial  coordinates  and 
then  quantized  in  brightness.  For  typical  images  (used  here),  256  *  256 

g 

image  samples  are  taken  and  the  brightness  is  quantized  to  2  levels. 

This  amounts  to  an  equivalent  data  rate  of  256  *  256  *  8  ■  5.24  x  10^  bits 
per  image.  Based  on  the  current  telemetry  channel  capacity  of  the  order 
of  16.2  kbits/second  [  10  ) ,  it  would  take  about  32  seconds  to  transmit 
5.24  *  10 **  bits.  For  many  other  data  sources,  such  as  satellite  and  bio¬ 
medical  image  sources,  the  data  rate  are  even  higher.  The  large  quantities 
of  image  data  generated  for  real-time  digital  processing  make  it  mandatory 
to  consider  data  compression  techniques.  Data  compression  refers  to 
representing  an  image  by  as  small  a  number  of  bits  as  possible  at  an 
acceptable  fidelity  level. 

Most  image  data  are  highly  redundant,  so  that  compression  can  be 
achieved  by  removing  this  redundancy  before  transmission  over  a  communi¬ 
cation  channel.  Data  compression  schemes  for  two- and  three-dimensional 
Images  have  been  studied  by  many  researchers.  One  particular  class  of 
methods  is  called  statistical  data  compression  where  a  class  of  images 
is  characterized  by  a  stochastic  process  which  is  described  in  terms  of  a 
mean  and  covariance  function  represented  by  a  mathematical  model. 
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Franks  [  12  1  has  suggested  a  model  for  spectral  density  of 
two-dimensional  video  signals  from  which  a  separable  covariance  function 
can  be  developed.  Nahi  et  al  [  48,49  ]  and  Habib  1  [  20  1  and  others 

have  used  this  separable  function  model  for  processing  images.  This 
model,  although  mathematically  convenient,  is  not  very  representative 
for  many  images.  Jain  [  31  1  has  proposed  an  approach  where  the  image 
models  are  obtained  by  a  finite  difference  approximation  of  second  order 
partial  differential  equations  (PDFs).  This  approach  offers  a  rich 
mathematical  basis  for  formulating  and  solving  many  important  problems 
in  image  processing. 

In  this  research,  there  are  several  objectives  to  be  achieved.  The 
primary  objective  is  to  further  investigate  the  two-dimensional  image 
models  developed  by  Jain  [  31  ]  for  intraframe  image  data  compression. 

The  use  of  image  models,  i.o.,  causal,  nemlrausal  and  noncausal,  shows 
the  connection  between  predictive,  hybrid  and  transform  coding  schemes. 
From  the  choice  of  an  appropriate  image  model,  one  would  be  able  to  deter¬ 
mine  an  efficient  algorithm  for  data  compression.  The  second  objective 
is  to  extend  these  modeling  concepts  to  images  which  nay  be  viewed  as 
a  composition  of  a  stationary  field  and  a  deterministic  field  (e.g., 
edges  which  may  vary  with  different  images).  Therefore,  a  coding  technique 
can  be  designed  to  the  sharp  brightness  changes  and  also  to  preserve  the 
nice  property  of  conventional  transform  coding.  This  technique  is  called 
feature  transform  coding. 
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There  are  several  important  data  compression  techniques  or  concepts 
relevant  to  this  research.  They  are  briefly  stated  as  follows: 

1.2.1  Differential  Pulse  Code  Modulation  (DPCM) 

The  basic  principle  of  a  DPCM  system  is  to  quantize  and  encode  the 
changes  between  successive  samples,  rather  than  the  instantaneous  sample 
values.  It  was  first  proposed  by  Cutler  [  9  ]  in  1952  that  a  reduction 
in  quantization  noise  for  the  same  bit  rate  could  be  obtained  by  taking 
advantage  of  the  redundancy  or  correlation  between  samples.  Such  cor¬ 
relation  exists  in  practically  all  signals.  It  can  be  shown  that  the 
variance  of  the  difference  signals  is  considerably  smaller  than  the 
variance  of  the  original  signals  unless  there  is  no  correlation  between 
successive  samples,  Harrison  [  23  ]  showed  the  DPCM  concept  through  the 
use  of  predict  ion  theory.  hater  work  published  by  O'Neal  (  50  ]  showed 
a  design  procedure  for  DPCM  and  its  application  to  television  signals. 
His  results  gave  some  verification  of  modeling  video  signals  by  a  first 
order  Markov  process.  Other  work  related  to  the  design  of  DPCM  systems 
may  be  found  in  [10,45,68].  In  our  work  we  will  use  DPCM  schemes  in 
conjunction  with  transform  coding  and  show  their  origin  in  some  of  the 
image  models. 

1.2.2  Transform  Coding 

An  ideal  transform  coding  system  consists  of  a  u.iitary,  linear 
mapping  of  a  block  of  correlated  samples  into  a  set  of  statistically 
independent  coefficients.  Data  compression  is  then  accomplished  by 
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sorting  and  quantizing  the  transform  coefficients  according  to  their 
information  capacity. 

Andrews  and  Pratt  [  5  ]  introduced  the  concept  of  transform  coding 

of  images  via  Fourier  transform  in  1968,  which  led  to  the  investigation 
of  utilizing  various  unitary  transforms  for  image  processing  [1,3,4,21,26,55,70] 
Among  the  candidates  for  transform  coding,  the  Karhunen-Loeve  transform 
provides  minimum  mean  square  error  for  any  fixed  rate.  Common  image 
transforms  are  discrete  Fourier,  Hadamard,  Cosine,  Sine,  Slant.,  etc. 

(refer  to  Chapter  4). 

In  our  work  we  will  show  that  a  particular  class  of  image  models, 
namely  noncausal  models,  lead  naturally  to  transform  coding  algorithms. 

It  will  be  shown  that  the  choice  of  transform  and  the  coder  design 
parameters  are  determined  via  the  model  parameters. 

1.2.3  Hybrid  (Transform/DPCM)  Coding 

Hybrid  coding  is  a  technique  which  combines  the  advantages  of  trans¬ 
form  coding  and  DPCM  coding,  which  was  first  proposed  by  Habibi  [  17  ] 
for  two-dimensional  images  modeled  by  a  separable  covariance  function. 
Subsequently,  Roese  et  al  (  58  ]  extended  this  technique  for  interframe 
coding  of  images.  Recursive  filtering  of  two-  and  three-dimensional  images 
based  on  hybrid  coding  concept  was  proposed  by  Jain  and  Angel  [  33  ] • 

Other  studies  of  hybrid  coding  schemes  can  be  found  in  [  18  ]• 

In  our  work,  semicausal  image  models  lead  naturally  to  the  hybrid 
coding  algorithms  and  the  adaptive  and  nonadaptive  schemes  of  implementing 
those  algorithms  are  presented  in  detail. 
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1.2.4  Contour  Coding  and  Edge  Extraction 


Since  a  significant  amount  of  visual  information  in  an  image  lies 
in  the  edges,  it  is  reasonable  to  consider  coding  schemes  where  the 
edges  are  treated  more  adequately  so  that  the  processed  image  would  appear 
acceptable  to  a  human  observer. 

In  1958  Schreiber  and  Graham  (15,62-65  ]  introduced  an  edge  coding 
technique  called  synthetic  highs,  in  which  a  scanned  video  signal  is 
segmented  into  its  high  and  low  frequency  components  and  they  are  trans¬ 
mitted  separately.  The  high  frequency  component  essentially  contains  the 
edge  information.  At  the  receiver,  the  high  frequency  component  is 
synthesized  from  the  addresses  of  the  image  edges  and  combined  with  low 
frequency  component  to  gener.-.ca  the  original  video  signal.  Yan  and 
Sakrison  [  74  ]  proposed  a  similar  scheme  where  an  image  is  represented 
by  two  sources,  texture  and  edge.  Image  quality  can  be  improved  by 
simply  preprocessing  with  an  appropriate  method  before  transmission  and 
postprocessing  wi th  its  inverse  after  reception. 

Use  of  image  models  for  segmenting  high  and  low  frequency  compo¬ 
nents  will  he  shown  in  this  research  where  several  simulation  results 
are  illustrated  to  show  the  effectiveness  of  the  new  process. 

1.3  Organization 

Having  outlined  the  background  material  and  the  importance  of  data 
compression,  we  can  now  describe  the  organization  of  the  thesis. 

Chapter  two  provides  a  general  discussion  of  image  models  and  their  pro¬ 
perties.  The  efficiency  of  these  models  is  measured  by  comparing  their 
data  compression  ability  in  terms  of  basis  restriction  errors  (for  defi¬ 
nition  see  Chapter  two).  Chapter  three  is  devoted  to  the  development  of 
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Hybrid  (Transform/DPCM)  coding  algorithms.  It  Is  shown  that  the  semi- 
causal  class  of  models  leads  to  such  algorithms.  Four  different  semi- 
causal  models  have  been  considered.  Adaptive  schemes  are  designed  by 
simply  adapting  the  model  parameters  to  cope  with  the  changes  in  image 
statistics.  The  effect  of  a  noisy  transmission  channel  on  the  Hybrid 
coding  system  is  also  studied.  It  is  shown  that  the  semlrausal  image 
models  may  be  used  to  develop  a  Hybrid  coding  system  for  compression  of 
noisy  images.  In  Chapter  four,  we  show  how  the  noncausal  image  models 
may  be  used  to  design  transform  coders.  Adaptive  and  nonadaptive  algo¬ 
rithms  are  considered.  In  Chapter  five,  we  extend  the  algorithms  of 
Chapter  four  to  represent  images  as  a  composition  of  two  sources.  One  of 
the  sources  is  a  stationary  component  and  the  other  is  a  nonstationary 
or  deterministic  component.  Examples  of  computer  simulation  of  all  the 
algorithms  on  image  data  are  given  and  comparisons  arc  made. 

Finally,  the  Appendices  contain  the  derivation  and  proofs  that 
relate  to  various  topics  throughout  the  thpsis.  In  particular,  Appendix 
A  provides  useful  tabulated  data  for  various  quantizers  and  serves  as  a 
good  reference  for  engineers  interested  in  designing  quantizers  for  image 
processing  applications. 
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CHAPTER  TWO 


DISCRETE  GAUSSIAN  RANDOM  FIELDS  AND  STOCHASTIC  IMAGE  MODELS 

2.1  Introduction 

In  discussing  image  representations,  we  think  of  characterizing 
images  as  two-dimensional  random  fields.  If  D  is  a  constant 
coefficient  difference  equation  operator,  we  consider  representations 
of  the  form 

Dlu^]  -  tti  (2.1-1) 

where  denotes  the  array  of  image  pixels,  (c^)  is  a  two  dimensional 

white  noise  process  or  a  moving  average  field. 

In  this  chapter,  we  consider  causal,  semicausal  and  noncausal 
representations  which  are  of  the  form  of  (2.1-1).  Those  representations, 
as  we  will  see  in  Chapters  3,  U  and  5,  will  aid  in  developing  useful 
algorithms  for  many  imago  coding  problems.  Here  we  will  study  data  com¬ 
pression  efficiency  of  these  models  measured  by  the  basis  restriction 
error. 


2.2  Causal,  Semicausal  and  Noncausal  Estimates 

Let  {u  }  represent  a  zero  mean,  stationary  two-dimensional  dis- 

»  »  J 

Crete  Gaussian  random  field  whose  covariance  function  is  defined  as 


E  Ui,JUi+k,J+l  "  *  (2.2-1) 

Let  u  denote  an  estimate  of  the  random  variable  u  .  We 

*  •  J  *  *  J 

consider  the  following  three  types  of  estimates. 


2.2.1  Causal  Estimate 

Suppose  the  elements  of  the  random  field  (u  ,}  are  arranged  in  any 

*  *  J 
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desired,  one-dimensional,  ordered  sequence.  Then  u 


is  a  causal 


estimate  of  u  if  It  depends  only  on  the  elements  that  occur  before 
1  •  J 

the  element  u  . .  A  common  example  occurs  when  the  image  is  scanned  column 
1 » J 

by  column  and  u  is  a  linear  estimate  based  on  all  the  elements 
*  *  J 

scanned  before  arriving  at  (i,j),  l.e.. 


“i.J  "  l  a(m*n)ui-m,J-n  (2.2. 1-1) 

m.nrS 

where  S  -  {m,n:  n  >  0,  V  ra)  |J  {m,n:  n  •  0,  m  >  0).  Fig.  2.1(a)  shows 
the  set  S  for  causal  estimate  at  (i,J). 


2.2.2  Semlcausal  Estimate 


If  the  estimate  u.  is  causal  in  one  of  the  coordinates  and 

1  •  J 

noncausal  in  the  other,  it  is  called  a  semlcausal  estimate.  For 


example,  a  linear  semlcausal  estimate  which  is  causal  in  "J"  and 
noncausal  in  "1”  would  he  of  the  form 


u.  .  -  }‘  a(m,n)u  (2. 2. 2-1) 

m.ncS 

where  S  -  {m,n:  n  >  0,  V  m)  I)  {ra,n:  n  »  0,  V  m  4  0)  and  is  shown  in 
FiR.  2.1(b). 


2.2.3  Noncausal  Estimate 

The  quantity  u  is  a  linear  noncausal  estimate  of  u  ,  if  it  can 
1 1  j  l » J 

be  written  as  a  linear  function  of  possibly  all  the  variables  in  the 
random  field,  except  u  .  itself.  For  example,  a  linear  noncausal 

*  •  J 

estimate  would  he  of  the  type 


(2. 2. 3-1) 


u.  .  -  l  a( m,n)u 

i-m,j-n 

m,ncS 

where  S  ■  {m,n:  V(m,n)  i  (i,j)}  and  Is  shown  In  Fig.  2.1(c).  Note  that 

u  contains  terras  from  all  the  four  quadrants  about  the  point  (i,j). 

1 » J 


2.3  Stochastic  Representations  of  Gaussian  Random  Field 

Let  u  be  an  arbitrary  estimate  of  u  . ,  then  a  stochastic 
i  •  J  * » J 

representation  of  the  Caussian  random  field  {u.  . )  is  defined  as 

*  »J 


"i.j  "  "i.j  +Ci.) 


(2.3-1) 


where  (c.  }  is  another  random  field. 

*  •  J 

There  are  two  types  of  represrntat ions  which  are  of  interest  to  be 
considered  here.  These  are  as  follows: 

1)  Minimum  variance  representation. 

2)  White  noise  driven  representation. 


2.3.1  Minimum  Variance  Representation 


A  minimum  variance  estimate  is  one  which  minimizes  the  mean  square 


'i.j  '  E<ui.j  -  "i./ 


(2. 3. 1-1) 


at  each  (i.j).  For  Gaussian  random  field,  the  minimum  variance  estimate 
would  be  linear.  The  coefficients  a(m,n)  In  (2. 2. 1-1),  (2. 2. 2-1)  or 
(2. 2. 3-1)  could  be  determined  from  the  covariance  function  r(V,t)  and 


the  orthogonality  relation 
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2.4  Spectral  Density  Functions  of  the  Gaussian  Random  Field  Represented 
by  Linear  Models 

Consider  the  stationary  Gaussian  random  field  {u  .}  defined  on  the 

* 

infinite  plane,  -*»  <  i,j  <  ">.  The  two-dimensional  e-transform  of  the 
covariance  function,  called  the  covariance  generating  function,  is 
defined  as 

*  Two-dimensional  e-transform  is  defined  as 

qp  ao 

r(z.tt7)  -  l  l  fd.j)*"1*"^. 

1  i—  j— 
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Su(z1,JE2)  “  I  l  r<k,£)r“kr2k 

Substituting  and  e^‘J;  for  r.j  and  z.,,  respectively,  we  have 


(2.4-1) 


S  (<u,  ,U)„)  ■  S  (z,  ,z_) 
u  1  2  u  1  2 


* ^  L  »  *2  "  • 


(2.4-2) 


This  becomes  the  spectral  density  function  (SDF)  of  (u.  }.  For  simpli- 

i  •  j 

city,  we  refer  to  Sy(z^  z^)  as  the  spectral  density  function  of  the 

Caussian  random  field  {u.  .)  at  various  places  In  the  following  chapters 

*  I  J 

without  further  not i f icat i on . 

The  discrete  operator  of  (2.1-1)  can  be  written  as 


DfZj  Zj)  -  1  -  l  a(m,n) z^r 
o.nr.S 


(2.4-3) 


If  S^fZj.Zj)  is  the  spectral  density  function  of  {c^  ^) ,  then  the 

spectral  density  function  of  fu  is  given  by 

1  *  J 


S  (z  ,z.) 

su(*i**2^  ’  rj  _i  • 

DfZj.Zj)  D  ( z x  ,z2  ) 


(2.4-4) 


The  covariance  function  of  (u  .}  is  obtained  simply  by  the  Inverse 

*  *  J 


Fourier  transform  of  S  (z,  z„) 

u'  1,  2 

2tt  2  it 


R(k,f.)  -  — *2  |  |  Su(«j,Cl>2)  exp(ju1k+J<»)28.)da>1d<iJ2  (2.4-5) 


0  0 


0 
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2.5  Image  Representations 


Now  we  consider  different  types  of  random  field  representations  as 
candidates  for  image  modeling. 


2.5.1  Causal  Representations 


The  general  form  of  causal  representations  of  an  image  field  is 
given  by 


Ui.J  "  Ui.J  +  Ci,J 


(2.3-1) 


where  u  is  the  causal  estimate  of  u.  ,  defined  in  (2. 2. 1-1). 

MJ  1  *  J 

As  an  example,  one  causal  representation  is  of  the  type 


*1,1  '  *lUi-l,l  +  Vl.j-l  +  Vl-l.j-1  +  C1.J  (2  5.X-1) 

whose  spatial  structure  is  depicted  in  Table  2.1  and  it  is  obvious  that 
u  is  causally  related  with  three  points  which  occurred  at  the  left 
upper  quadrant.  It  has  been  shown  that  this  equation  and  other  similar 
causal  representations  can  be  obtained  by  discrete  approximation  of 
hyperbolic  partial  dlfferentail  equations  [  34]. 

Using  the  spatial  structure  and  measured  covariances  of  a  given 
image  field,  one  may  identify  the  coefficients  a^,  a^,  a ^  and  the  random 
field  {e.  .)  via  least  square  regression  techniques.  For  the  covariance 

MJ 

function 


o \  *  Ikl  It! 

r(k,t)  -  EuijUl+k>J+l  -  Pl  p2 


(2. 5. 1-2) 


where  p^  and  p2  are  the  vertical  and  horlsontal  correlations  respectively. 
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The  model  of  (2. 5. 1-1)  is  its  exact  realization  when 


al  "  P1  *  a2  “  p2  »  a3  *  “P1P2 


(2. 5. 1-3) 


and  (c.  .)  is  a  white  noise  field  with  zero  mean  and  covariance  function 
*  •  J 

E  '  B\.0{1.0  •  »2  ‘  (I-PJHI-P2)-  (2- 5.1-4) 


The  corresponding  discrete  operator  and  spectral  density  function  of 
(2.5. 1-1)  subject  to  (2. 5.1-3)  are  given  as  follows: 


DCz^Zj) 


Sc(z1,72) 


(l-pl*11) (l-p2z21^ 
e2  -  (i-p2)d-p2) 


(2. 5. l-5a) 


(2.5.1-5b) 


Su(zrz2) 


(1-PjZj) (1-pjZj1) (l-p2*2) (l-pj*'1) 


(2.5.1-5c) 


In  image  processing  literature,  this  model  has  been  used  for 
differential  PCM  and  filtering  of  images  [34]. 


2.5.2  Semicausal  Representations 

A  semicausal  representation  is  causal  in  one  of  the  image  coordinates 
(say,  j)  and  is  noncausal  in  the  other  (say,  1).  For  example. 


“i.i  '  “<ui-i.)  +  “i+i.j'  *  Y“i.i-i  +  ci.j 


(2. 5. 2-1) 


where  Ial  <Js  ,  lYi  <  1  and  j  2a  +  Yl  <  1. 


(2.5.2-la) 


This  is  a  semicausal  representation  of  the  Gaussian  randoai  field. 
Conditions  of  (2.5.2-la)  insure  the  stability  of  (2. 5. 2-1).  The 
covariance  function  generated  by  this  equation  depends  on  the  covariance 


« 


function  chosen  for  the  random  field 
is  a  white  noise  field,  i.e.. 


). 


For  example,  when  (c 


i.J 


} 


rc(k,i)  -  Eri>jc1+vCtj  +  )l  "  8  6k,06t,0* 


(2. 5. 2-2) 


The  spectral  density  function  of  {u.  }  is  then  given  by 

i  *  J 

.2 


Su(zl’z2) 


0 


-1  -1  "I  * 

(1-QZj-aZj  -yzjKl-aZj  -az^yzj  ) 


(2. 5. 2-3) 


This  white  noise  driven  representation  will  be  called  the  white  noise 
driven  SCI  model. 

As  a  generalization,  consider  the  case  when  {e.  .)  has  its  covariance 

*  t  J 

function  given  by 


rc(k,fc)  -  Eet  jc1+ktj+t  "  8  6i,0(5k,0"otl6k-l,0"al6k+l,0)'  (2-5-2_A> 


The  corresponding  SDFs  are  given  by 


Sc^Zl,Z2^  "  P2l1_<l1Cz1+zi1)  1 


Su<*l>*2)  ' 


(l-az^-ar^-yzj)  (l-az^-az^yXj1) 


(2.5.2-5n) 

(2.5.2-5b) 


For  »  0,  (2. 5. 2-3)  results,  we  have  white  noise  driven  model. 

For  tij  •  o,  the  quantity  (see  (2. 5. 2-1)] 

ui,j  "  °(Vl,j  +  ui+l,j)  ♦  (2. 5. 2-6) 

becomes  the  best  linear  mean  square  estimate  of  u  given  k,i  belonging 

+  •  J 

to 

S  -  (k,t:  l  <  J,  V  k)  U  {k,t :  l  -  J,  V  k  i  i) . 

This  follows  by  observing  that  the  orthogonality  relation 
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I 


■  "t.J> 


■  Eui.)cl.)  ■  “• v  1>J 


(2. 5. 2-7) 


is  satisfied  when  o.  ■  a. 

Then  (2. 3.2-1)  becomes  a  minimum  variance  SCI  representation,  since 

*  2 

E(u  -  u.  )  is  minimum  over  the  set  of  lattice  points  of  S.  The 
spectral  density  function  of  the  minimum  variance  semicausal  representa¬ 
tion  is  obtained  by  setting  ■  a  in  (2. 5. 2-5). 

It  has  been  shown  by  Jain  [  28  )  that  the  separable  covariance 
function  [see  (2. 5. 1-2)]  has  a  semicausal  minimum  variance  realization 
given  by  • 


Ui,j  ’  al(ui-l,j+lii+l,J)  ~  P2Cll(“i-l,j-l+ui+l,j-l)  +P2Ui,J-lCi,J 


where  (c  }  is  chosen  to  have  its  covariance  given  ns 
*■  t  J 


(i+pJ) 


£.,0'  “luk,-l  ~rk,l 


"k,0' 


(2.5. 2-8) 


(l-p?)(l-p*) 

r^(k,£)  -  - - — ^ — —  6,  n(~a,S,  .  -  a,6,  ,  +  6,  n>  (2. 5. 2-9) 


and 


Ql  “  .  2 

1+pl 


with  corresponding  SDK 


(l-p?)(l-p*)  . 

S  (z  z  ) - - -  (l-o  z  -a  z  ) 

1  (1+pJ)  1111 


(2.5.2-10) 


This  will  be  called  the  minimum  variance  SC2  representation. 
Interestingly,  the  spectral  density  function  S^z^.z^)  of  this  minimum 
variance  SC2  representation  and  that  of  the  white  noise  driven  Cl  repre¬ 
sentation  are  identical.  For  more  general  minimum  variance  semicausal 
models,  see  [  34  ] . 
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The  semicAU9al  representations  will  be  shown  to  lead  some  useful 
hybrid  algorithms  for  image  coding.  A  hybrid  algorithm  is  a  recursive 
algorithm  obtained  after  transforming  each  row  (or  column)  of  the  image 
by  a  unitary  transformation.  Note  that  above  mentioned  semicausal 
representations  can  also  be  obtained  by  discrete  difference  approximation 
of  parabolic  partial  differential  equations  I  34  ] . 


2.5.3  Noncausal  Representations 

The  noncausal  representations  are  boundary  value  problems  in  each 
coordinate  and  have  correspondence  with  the  elliptic  class  of  partial 
differential  equations  (PDEs)  [  34  ] .  For  example,  the  noncausal 
representation  (NCI) 


ui,j  "  a(ui-i,j  *  Vi.j  +  ui,j-i  +  Ui,j+1)  +  ci,j  ’  1“) <!< 

(2. 5. 3-1) 

is  a  discrete  spproximat ion  of  a  Poisson  PDE . 

In  general,  (e  .  ^  may  be  chosen  9uch  that 
1  *  J 

rc<k,l)  •Ec1>jti+k,.)-fk"e  l5k,o<l,o"al<{k,-l*Sk,l>5t.O 

•  |ai'  ‘  •>  <2-5-:>-2) 

The  spectral  density  functions  for  {c.  ,}  and  fu  .}  are  then  obtained  as 

*»J  i»J 

(2. 5. 3-3) 

*  A 

Sc^*l ,b2> 


Sf (*i »*2^  “  “  a1(E1+*i1+*2+r21^ 


•u<*r,2> - 


[1  -  a(z1+s11+t24-x“1)]2 


(2. 5. 3-4) 
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There  are  three  types  of  NCI  models  to  be  considered  here. 

a  -  0.  {c  .)  becomes  a  white  noise  process  and  (2. 5. 3-1)  is  a 
i  i  *  j 

white  noise  driven  NCI  representation. 


-  a,  we  get  the  minimum  variance  NCI  rep 


i.  Note  that  In 


the  spectral  density  function  a  pole-zero  cancellation  occurs 

and  the  S  becomes 
u 


Su(zl’*2) 


[1  -  a(’1+z”1+z2+z“1) ] 


(2.5. 3-5) 


c.  For  cij  i  a,  (2. 5. 3-1)  may  be  considered  as  a  noncausal  autoregressive 
moving  average  (ARMA)  model,  denotes  a  NCI  ARMA  model.  According  to 
Jain  [  341,  model  identification  experiments  have  shown  for  several 
images  (for  example,  Cirl,  Moon  and  Earth)  that  this  model  gives  a 
better  fit  to  the  real  data  than  the  previous  two  models. 

Another  type  of  noncausal  representation  can  be  derived  from  a 
bihanaonic  PDE  (34’,  and  gives  a  13-neighbor  model 


ui,j  "  a(Vl,j+Vl,J+Ui,J-l+Ui.J+l)‘aa<Ui-l,j-l  +  Ui-l,j-Fl 

e 

+  u  i+l,J-l  +  ui+l,j+l)  *T(ui-2,)+ui+2(J+ui,j-2+ui,J«)  +  (l44a2} 

(2. 5. 3-6) 


where  a  ■ 


(l+4a2) 


If  covariances  of  {e  .}  are  taken  as 

I  •  J 


1 
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rc(k,£) 


Ec  c 

“-i,jcl+k,j +1 


B2  < 


-al 

,  if 

It-  ±1, 

k-0  or  it-0. 

alal 

.  if 

It-  ±1. 

k-  ±1 

alal 

2 

.  if 

*-±2  , 

k-0  or  k-  i2 

1 

.  if 

k-£.-0 

0  ,  otherwise 


(2. 5. 3-7) 


where 


iv  <  5  »  iaii  <  a  ; 


then,  for 


a^  -  a  ,  -  a 


Eqn.  (2. 5. 3-6)  becomes  a  minimum  variance  NC2  representation.  Spectral 
density  functions  for  this  mode  are  given  by 


Se(V*2>  "  B2[l  -  ajtzj+z^+zj+z”1)]2 


Su(al-*2)  “ 


Sr  (Zj  ,7-2) 


(1  -  aCz^z^+z^+z-j1) 


(2. 5. 3-8) 

(2. 5. 3-9) 


The  noncausal  representations  discussed  above  often  realize  nearly 


isotropic  S^Zj  z2), 


Not  all  noncausal  models  give  near-isotropic  random  fields.  For 
example,  consider  the  noncausal  model  (NC3) 


Ui,j  "  “l^i-l.J  +  Ui+l,J)  +  a2(ui,j-l  +  “i.j+l* 


(2.5.3-10) 


“  aia2(ui+l,j+l  +  Ui+l,J-l  +  Vl.J+l  +  Ui-l,j-l)  +  Ci,J 


If  the  covariances  of  (g.  ,)  are  given  by 


1 


if  k-e-0 


rc(k.t)  -  fcljEitkj+t  -  62  1 


-°1  • 

if  k*>  ±1,  i-0 

"a2  * 

if  k-0,  i-  ±1 

V2  ' 

if  k-  ±1,  i- 

0  , 

otherwise 

(2.5.3-11) 


(l-p?)(l-p*) 

where  B  -  - = - r~ 

(l+ppCl+pp 


a  ~ 

1 


(1+Pj) 


2  a+Pj) 


SDFs  of  (2.5.3-10)  are  given  by 


SC(Z1,Z2)  " 


B 


(1  -  OjUj+z,  *)  ]  [1  -  a1(z2+z21)] 


(2.5.3-12) 


Su(zl-z2)  ‘ 


(l-pj)(l-p2) 


(1-PjZj)  (1-PjZj1)  (l-p^Ml-p^1) 


(2.5.3-13) 


Notice  that  the  form  of  S  (z,z„)  is  the  same  as  the  SDK  of  causal  DPCM 

u  1  2 

Cl  model  and  semlcausal  SC2  model,  so  it  generates  the  separable  covariance 
model  which  is  non-isotropic. 


2.6  Basis  Restriction  Data  Compression  Efficiency  of  Image  Models 

The  performance  of  the  foregoing  image  models  in  image  data  com¬ 
pression  can  be  analyzed  by  studying  their  basis  restriction  efficiency. 

The  scheme  is  shown  in  Fig.  2.2.  First,  the  original  image  (u  }  Is 

i  •  J 

transformed  by  a  unitary  matrix  A  (see  Chapter  4,  the  definition  of 
transformation).  Next,  a  two-dimensional  sample  selection  function 
{w(i,J)}  (also  called  zonal  filter),  which  at  each  (i,J)  takes  on  the 
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value  of  zero  or  one,  is  applied  to  the  transformed  image  (v.  .  Only 

*  •  J 

those  samples  for  which  (w(i,J)"l}  are  stored,  zeros  placed  at  other 

locations.  This  zonal  filtered  image  (v,  .)  is  inverse  transformed  by 

1 » J 

A  to  give  { u  .) . 

^  •  J 

The  basis  restriction  mean  square  error  (B.R.M.S.E.)  is  defined  as 


B.R.M.S.E. 


'  \\  E  \  < 
'l]*  (vi.j  -  'i.j  )2/\  \  vij 


(2.6-1) 


where  we  have  used  the  fact  that  A  is  a  unitary  transform.  Associated 
with  this  error  is  the  sample  reduction  ratio  (S.R.R.) 

Number  of  image  samples 

S.R.R.  - - .  (2.6-2) 

Number  of  samples  selected  by  zonal  filter 

This  ratio  is  a  measure  of  data  compression. 


2.6.1  Zonal  Filter  Design 

The  zonal  filter  w(i,j)  is  designed  as  follows. 

1.  Let  t  be  the  number  of  samples  to  be  selected  for  transmission  for 

n 

a  given  S.R.R.  Clearly  t  “  NM/S.R.R.  for  an  M  *  N  image. 

2 

2.  Find  the  variances  o  of  the  transformed  image  samples  v  .  If 

m,n  **  r  m,n 

2 

A  is  the  DFT,  and  N  is  large,  an  estimate  of  o  is 
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2  I 
a  -  { 

m,n 


2 

0  H-o,N-n 


0<  !»l  i  I  .  0<  |n|  <  H 

othervise 


3. 


will 


2 

Arrange  a  in  a  decreasing  order.  Then  w(m,n)  -  1  for  the  first 
m  v  n 

t  addresses  (m,n)  and  w(m,n)  “  0  othervise. 
n 

2 

Since  w(m,n)  depend  on  0  or  the  SDF,  the  zonal  filter  functions 

m,n 

be  different  for  different  image  models. 


2.6.2  Experimental  Results 


For  easy  reference,  all  the  image  representations  (models)  presented 
in  Section  2 . S  are  summarized  in  Table  2.1.  It  is  seen  that  various 
partial  differential  equation  (PDF)  operators  are  also  listed.  Table  2.2 
shov.  model  parameters  used  for  three  different  images.  These  parameters 
were  found  by  least  square  matching  of  the  16  *  16  image  covariances  with 
the  model  covariances. 

Simple  basis  restriction  error  experiments  have  been  performed  to 
compare  the  various  imago  models  for  the  Girl  and  Farth  images.  We  have 
used  the  discrete  Fourier  transform  to  be  the  unitary  transform  operator. 
Since  the  image  models  differ  in  their  SDFs,  in  the  transform  domain. 


the  variances  of  the  samples  v  will  have  different  distributions. 

*  »J 

Figs.  2.3  -  2.5  show  the  zonal  filters  at  various  sample  reduction  ratios 
for  causal,  semlcausal  and  noncausal  models.  Table  2.3  lists  the  results 
of  basis  restriction  mean  square  error  versus  image  block  size  for  the 
Girl  linage.  Table  2.5  lists  the  corresponding  results  for  the  F.arth 
image.  Figs.  2.6  -  2.7  show  the  plots  of  basis  restriction  errors  versus 
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sample  reduction  ratios.  In  general,  basis  restriction  error  of  the 
noncauaal  model  is  lower  for  most  sample  reduction  ratios  (indicating 
its  superior  fitting  of  the  image*  data  than  the  other  two  types  of  models. 
For  very  large  sample  reduction  ratios,  since  few  transform  domain 
image  pixels  are  transmitted,  all  the  models  select  the  same  pixels 
(see  Figs.  2.6  -  2.7)). 

The  above  conclusion  could  also  be  drawn  from  Figs.  2.8  and  2.9, 
which  show  the  Basis  Restriction  Data  Compression  of  Clrl  and  F.arth  images 
using  a  16  *  16  image  block  size.  Clearly,  noncausal  models  give  the 
best  results  followed  by  semicausal  models. 

2.7  Summary  and  Conclusions 

In  sunsnary,  several  image  representations  have  been  considered. 

They  originate  from  finite  difference  aoproximat ions  of  certain 
PDEs.  These  models  have  been  compared  on  the  basis  of  their  data  com¬ 
pression  efficiency  based  on  a  basis  restriction  error  criterion. 

The  following  conclusions  are  made. 

1.  The  models  Cl,  SC2  and  NC3  all  generate  random  fields  whose 
covariance  functions  are  the  same  as  the  separable  model  [see 
(2. 5. 1-2)]. 

2.  For  cij  »  a,  the  spectral  density  of  NC2  is  the  same  as  that  of  NCI 
for  Qj  -  0.  However,  NC2  is  a  minimum  variance  noncausal  represen- 
tion,  whereas  NCI  is  a  white  noise  driven  model. 

3.  In  a  data  compression  application,  the  noncausa!  models  give  the 
best  performance.  A  similar  conclusion  has  been  made  In  image 
filtering  applications  [34]. 
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Figs.  2.3  -  2.5  could  be  Interpreted  as  the  spectral  density 
function  Indicators.  Different  types  of  image  models  would  achieve 
different  shapes  of  spectral  density  functions.  Figs.  2.3  -  2.5 
show  the  contours  of  approximately  equal  spectrum  values. 

Higher  order  models  may  be  needed  to  represent  image  SDFs  more 
accurately.  Identification  of  such  models  remains  an  open  problem. 


23 


3ot>:oc  ^a.'oom::* ' *"ioc: **  o.o 
3:0:0000300  .  mime;.?'*. ..  -to 

3i?3oc>cooioen.  *.::.::*  .*'•»•.  :•  jo 
3U3*.  .u3ut.D..l4*l.*..  ..'•*•*  -  .  .  C 
oootvco:  0:33: m: m: noaio.  *.000 
o-  .'coeoceacomii  in: .  joro.-'  v.'oo 

3.1';  .t  ioocjicx;:: l _  : .nut  cn: 

ooooooo,.*c3..ii:im . .  ■  -'cot*: .  uoro 
oo3joc33joonu:i. : . . : n;?o'  . coo 
OCOC03003 iiuuui n: :  occoocooo 
ooooocoooiuiiuiimiiiocc'ivooo 
00030031  lull  i.ll  IV  ill  t  UlvOuO'JO 

ococi ;  m  i  m  h  u  i  i  1 1  u  i ;  1 1 1  iooo 
mmiuimiumimimiuu 
mmuiiaumnimiunuu 
m  i  i  1 1 1 : 1 1 1 1  i  i  u ;  1 1 1 :  u : :  n :  n  i 
nmmimiimnimmuun 
lmuuum  minim  miiim 
linn  mini  mi  mini  mi  mi 
1 1 1  u  1 1  m  l : :  1 1 1 » u  u  u  l  i  i :  l  m  l 
eeioimiiiii  mi  liiiiiimi  ioco 
OOOo.'COU  111111111  mill  1900  0000 
00093C0C3 111 1 1 1 111 : 11 11 1 0C00  3300 

oossooooooimmmmccooccooo 

033G3C0C3J11 1 11 1 1 i 1 11 10C33030300 
00303C0C0OD . I  11111  111 33000CC COCO 
00000030030311:111: UOCOCOOCCOOO 
003  030200  3301  nil l:  110000333 33  30 
00003000330U111: 11 l 1103033 CO 03 30 
0300 33000300 Oil il i 1 130330333 fCPC 

oooojcoco.'ccc:  1  nil  looocooooocoo 
003300330333c:: i: : ncaocoococooo 

(a)  3:1  Saaplc  Reduction 

O0C0300O3 3C000C 3 1 0330030  >'.3C<':30 
03 330000033 30C0C 10 30003 C30C 3 COC'9 
00 03 03 00C 3 330 300: C  330030*30  33C3C3 
aoooccoc-.' 33303301  333  303C 3  000  C300 

00000000303013.01  33 303  JC330CT  'Ov, 
000030 103 3 3300331 03 3003300COC333 
000C0C30J3C0J00310 J030303JOO  330 
00003CC33330’,  OC3:Dj33r3&33CO.COO 
OOOCPC3O33  3wJ0C01C33:JO  J.'Ct'CCCOO 
00O03  3C.  JD3C00011  l  3 0333 03003  tO 30 
0300CC30030C3001: lOOCGCuOOOdsODO 
0COC003003333U01 i . 300333 3033 1330 
030000000  33  3*3001  U  33003 OOCOCCOOC 
000300303  J3000nn:  J003O33C0C30C 
00003C03333001i:U : IOCjOwOCOCOOO 
30033033’: 11111  1 1  1 11111 103000000 

1 1  n  i  m  i :  1 1 1 1 1  n  1 1  u  n  1 1 1 1 1 1  u  i 
aoocooacoiiiniii  mill  iocoo  co  30 

OOU0033303330U  i  1  1  *.  1  JC  3000u3  3C0  J 
00300C330 330331  Hi  : 0333 03000 3300 
0000033033330031 11 OOPC 30PCCO C300 
0030rcC0033C*3(l01  .  .  :01330?30'OCC3l 
0S30033CC330O30 l 1 1 OOuO 3 3330C CO 00 
OO0C3C3CO33C03O1  *  1  0033  JCC-P30  COCO 
30030  00C  3  33U  >i  3C0 1 00  OC  3 1CCCC3  tOOO 
OUOOCOOOO 33 0033 31  0  33 333033  IOCOOO 
003000333  300003  31  33.**  3330  JtOC  I'OCO 
03  303COC3  3  3C30(  J133C303C303CC030 
0  300*3300?  3  33C  0331  33  OJC  3  0 '.003  UC30 
003030003  130033  JiCC  *3033  003. 3  MCO 
0O<'00i30C00t30331OCPCC3C3C'1orCCO 
C0030C3CO  33ut033  tOOOiOOOOCCO  CO 00 
(c)  8:1  Snpl<  Reduction 


3050'.C3<>-,  j<*.*  *  >’.l  1  >.  )''313C'j. '*.'.'*.  30 
5333 30COD9C . DC  1 1 1 1 CuOO3003CO  C3  33 
000030  30*330003011  lOOCOOCOOOO COCO 
3  C  0 . 3  :  3 . C  330131*1  lOO.OOOOCCw  . 
00000 w  3.105  JCoO*. i  a  1'JoO  3  3  jOCOOOO  > 

00*3? o*j  i.‘i.*33'. .  jo i :  i  jcvcocosoj:.  .. 

300*3.3330  J3C0301 1  i  JO  JOO  30*103  uO  53 
3303000*.  33303301 1 1303330  COCO  COCO 

oocoooouoocco'ji:  1 1 io. 3 jcoouo to  ;o 

OCCCOCOC333CCOmilC3CJOO\<COCCCO 
30  03  5  COOO  30  f.t*  <Jli  11  *  ..33000COt.CC 

30000J  ,’coo  ;co;m  u  iooooooootcoo 
scoot  so.  ioool : i :i : i loooooocstooo 
oooocoDoooamm  nm.j  joococrco 
ooooccooi  mm  :n  nmiioooccoo 
iiiiiuiiuiiiiiiiuniiiiuii:i 
iimiimiiiiimi-iiimniiii 
mi  nnn  mi  m  mini  nn  n  n 

ooocoooou  i  m  1 1 1 1 1 :  M  i  iicoctooo 
ooooooooo'.'-3mm  i:ui3oo3oo  coco 

OPOOCOOCOIODCil’  1  U1300C3&OCC300 

ococ  rocooocoo:  m  moo  jcoooocc-'.o 

0000 000333 30 001  111 .C303CC0OC COCO 

ococ-cooo  Viooce:  1  •  1 100c  3o:ococc?3 
ooococ  ocooocjoi : :i i oooooccoccsc; 

OOOOCOOcCiOOOuOli iOODODOOCOOCOOO 
OOOOOC3053COJOO: 1 1 300030000000 :o 
030 0333U33 3 CfOCll 1 3 3333 033C0 CO 30 
00  300C1O2  JUCCJOU  iOOOCJCOCOOCCCC 
003300*3033  3  JCOOII  l  3 00 C  -  00ft  0  COCO 
300 300CC 330  0  :  30 1 1 1 503C0C 3330 COC3 
0002  JC CC DODO i*3Cl  1 1 00DC93C jOOCODP 
(b)  1:1  Sample  Reduction 

0000300C  "COJ  JOCC  30-330030  3000  1  0*0 
00000  3300  ;3Ct’CC31  0C3303 COCCOl  OC? 
OCOVO-’diCOOCOCOl  30033 3CC000 3330 
00  3  3  ;  u  JOC  30000001 330  )03.  3300 tOUC 
00 JOCOOOvtCOOOOOlOOOOCOCOOCOCCCO 
3COO0  33t3.*JtCCOCi03O30Jte.*C0CO0C* 
30000COC30C3C. 303 1 C0CC33.3CC0 COCO 
030  DOC  33030000  301  00  3303  CvO  JO  30*10 
3000C  003  330  300C0 1 03  .0 JOJCCC CO CO 
OOOOOOOCCCCOOTOOl 30J30 lUOOOCCSlO 
0C003C  ;OC?OOCOOOiCO;?03t*COOO?5?C 
000 30 COCO C3 30003 1030 DOC COOOC CO 33 
0003c033no0ucc00:03uc0000coot030 

0000CC3C0C33COC01 30D03 300 JCO 0033 
OOOCOOOOO’CCO 3001 OOOC 33000000030 
OOOOCC30  330COOJ11 13*330  330C00CC:3 

oimmimimuimiiimmi 

OOOOCOOCOCOtt 00 U ICOOOOOOCODOCOC 
OOCOOOCOC 2 3CC0C 31 035300 3C000 0300 
OOOOOCOtDOOCtuOOl C350C JC0CC3C003 
000000300 0050.0 01 CCOOGJ 330001 3 3C 
00000C0C3C03  30C01 03000  300003  COCO 
03  0  300000 ;0000031 COOCO ) 00500 tC03 
000003CCCC 3300031 ODDCDJ 30000 OC 33 
33OO053 OClOCuOOOi  OOOC'OCCCOGO  tv  DO 
OOOOOOOCOCOC 500 31  COO 3C 003000  0333 
OOOOOCOCOCO.OCC.'lCOCOOSOOOCGCOtC 
OOOOCOOuCCOCOOCOl  90030  3 0C900  COCO 
OOOOOUCOO 3000 DO  3 13000090 COOOC33C 
OOOOOOOCGOOJOOOCIGOCOOOOOOPOOOOO 
0033C330C*'  3C00031  00390  JC 3000  *.*3 00 
OOOCCCOOC.IOOPDOoIOOOO&OPCOOOOCcO 
(d)  16:1  Sample  Reduction 


Figure  2.1  Cauaal  Model  (Cl)  Fourier  Tranaform  Domain  Model  Zonal 
Fattema  for  a  32  *  32  Image  Block  Slae  with  Varioua 
Sample  Reduction  Rat  toe. 
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Figure  7.4  SnlciuMl  Nodcl  (SCI)  Fourier  Tranafore  Domain  Model  Zonal 
Fatterna  fnr  a  12  *  12  Inage  Block  Site  with  Virloua  Sample 
Reduction  Ratine. 
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(d)  16:1  Sample  Reduction 


Figure  2.J  Ko»rau«al  Hr*del  (!*C1  or  HC2-A  or  Fourier  Tranafona 

Domain  Hod el  Zonal  Pattern*  Cor  a  32  *  32  Image  Stock  Site 
with  Varloum  Sample  Reduction  Ration. 
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Figure  2.9  Fourier  Transform  Model  Zonal  Sampled  Karth  Image  with 
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CHAPTER  THREE 


SEMICAUSAL  MODELS  AND  HYBRID  CODING 


3.1  Introduction 

In  this  chapter  we  consider  Image  coding  techniques  which  combine 
transform  and  predictive  coding  techniques.  Such  techniques  have  been 
called  Hybrid  coding  |  17  ]  or,  more  generally.  Hybrid  processing.  The 
basic  idea  behind  hybrid  processing  is  to  first  convert  a  two-dimensional 
discrete  random  field  (an  array  of  random  variables)  into  a  sequence  of 
one-dlmenslonal,  independent,  discrete  random  processes  via  a  unitary 
transformation.  Then  each  of  these  one-dimensional  processes  is  pro¬ 
cessed  independently  by  one-dimensional  techniques  such  as  recursive 
filtering  (for  image  restoration  aopl icat ions)  or  predictive 
coding  (such  as  DPCM,  for  image  t ransmission/storage  applications). 

In  this  chapter  we  show  the  senicausal  models  considered  earlier 
yield  hybrid  coding  algorithms.  Once  the  parameters  of  a  semi  causal 
model  are  identified,  the  entire  hybrid  coding  scheme  can  be  designed. 
Adaptive  hybrid  coding  schemes  can  the"  he  designed  by  simply  adapting 
the  semi  causal  model  to  changes  in  image  statistics.  Some  practical 
adaptive  hybrid  coding  methods  which  strike  a  balance  between  complexity 
and  performance  are  studied.  Application  of  these  models  to  coding  of 
noisy  Images  are  also  considered. 

3.2  Scmicausal  Image  Models  and  Hybrid  Coding 

Let  tUj  represent  a  r.ero  mean,  Causslan,  stationary  random  field 
whose  covariance  function  is  defined  as 

'»•*>  "  *  ».!• 
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Consider  an  N  *  M  segment  of  this  random  field  where  1  <  i  <  H  and 


1  <  J  £  M,  and  let  u^  denote  the  Jth  column  of  the  Image,  i.e.. 


“j  '  (ul.Ju2.j- 


(3.2-2) 


From  (3.2-1)  it  follows  that 

EujuJ+l  ’  Rf.  (3.2-3.) 

where  R,.  is  an  N  *  N  Toeplitz  matrix  of  elements 

R^(m,n)  ■  r(m-n,l)  1  <  m,n  £  N  .  (3.2-3b) 


3.2.1  Hybrid  Coding  of  Separable  Covariance  Image  Fields 

Consider  the  often  used  separable  covariance  model  for  images 


r(k.f) 


•w 


(3. 2. 1-1) 


where-  p  and  p.  are  the  vertical  and  horizontal  correlations, 
v  h 

respectively. 


This  gives 


(3.2.1-2.) 


where 


(p 


»-n. 


(3.2.1-2b) 


It  is  well  known  that  Ry  is  the  covariance  of  a  aero  mean,  first  order 
stationary  Markov  process  whose  one  step  correlation  is  Py.  Let  ▼  denote 
the  Karhunen  l.ocvr  (KL)  transform  of  this  process  (see  Chapter  4),  i.e.. 


36 


y  is  a  unitary  matrix  that  satisfies 


yryt-a  yvT-i 

V 


(3. 2. 1-3) 


where  A  ■  {  X  is  the  diagonal  matrix  of  eigenvalues  of  R^.  Define  the 


t  runs  format  ion 


A 

vj  -yuj 


(3. 2. 1-4) 


Using  (3.2.1-2a),  (3. 2. 1-3)  and  (3. 2. 1-4),  one  obtains 


BwjvJ*i"fRivT"o2phlA 


EVj(i)vj+fc(i+k)  -  a2p;/ Vk,0 


(3.2.1-5a) 


(3.2.1-5b) 


Eqn.  (3.2.1-5b)  implies  that  1)  the  elements  of  the  transformed  column 
Vj  are  mutually  uncorrelated  and  ii)  for  each  i,  the  sequence  (v^(i), 

1  <  J  <  H}  Is  a  first  order,  stationary  Markov  process.  Therefore,  we 
can  write 


v,(i)  •  P.v  .(1)  +  e  (i)  (3.2.1-6a) 

J  h  J-l  j 


where 


Ee^(t)  -  0  ECj  (i)c  j+t(i+k)  -  O^l-p*)*^  Q6k  (3.2.1-6b) 

and  (v^(i),  1  <  i  <  N)  are  the  elements  of  the  vector  .  (3. 2. 1-6) 

spells  the  so  called  Hybrid  coding  method  for  the  random  field  (u  . 

»  »  J 

First,  each  column  u^  is  unltarily  transformed  (Ideally  by  the  KL  trans¬ 
form)  to  obtain  and  for  each  sequence  (v^(i),  J  "  1,2,...)  for  1  <  1  <_  N, 
an  Independent  DPCM  channel  is  used.  This  overall  scheme  is  shown  in 
Fig.  3.1.  Fig.  3.2  shows  the  details  of  the  ith  DPCM  channel  and  the 
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reconstruction  filter  at  the  receiver.  Often,  a  fast  unitary  transform 
is  used  to  replace  the  KL  transform  to  make  the  scheme  more  practical 
(see  Chapter  4).  For  video  image  data  where  P^  “  Py  =  0.95,  the  Cosine 
transform  performs  very  close  to  the  Kl.  transform  defined  in  (3. 2. 1-3) 
11,31,34] 

More  generally,  if  any  two-dimensional  random  field  (u  ,}  can  be 

•  J 

transformed  via  a  unitary  transformation  to  a  sequence  of  independent 


one-dimensional  Markov  processes,  then  a  Hybrid  coding  scheme  for  trans¬ 


missions  of 


} 


can  be  designed. 


3.2.2  Semlcausal  Representations  as  Two  Source  Models  for  Finite  Random 
Fields 

Consider  the  given  image  to  be  an  N  *  M  segment  of  the  infinite, 
stationary  random  field.  Semlcausal  representations  we  have  considered 
in  the  previous  chapter  (see  (2. 5. 2-1)  and  2. 5. 2-9)]  when  written  in  an 
N  *  1  vector  notation  are  of  the  form 
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(3.2.2-2b) 


1 


*J  “  1 


I N-2  termsi  T 

®Iuq  j  i  o,  ...  j  ^  SCI 


(3.2.2-2c) 


«(b  0 . 0,b  ]  ,  for  SC2 

^  I N-2  terms, 


bi.j  "  uo,j  "  PhU0,J-l. 


3N , j  "  Vl-l.j  "  PhUN+l,j-l 


(3.2.2-2d) 


The  N  x  1  vector  b^  contains  only  two  nonzero  entries  which  depend  only 

on  the  boundary  variables  { un  J  and{u  .)  of  the  1th  image  column. 

U,J  n+l.J 

Note  from  (3.2-2)  that  the  Jth  image  column,  contains  only  N  elements. 
The  covariance  of  tj  is 


“c  ■  EtjS  ' 


B2lS  ,  for  SCI 

J  I*1 


B  Q6  v  *  for  SC2 
j 


(3.2.2-2e) 


where  02  -  (1  -  p2)(l  -  p2)  /  (1  +  p2)  for  the  SC2  model.  (3.2.2-2f) 

v  h  '  v 


F.qn.  (3.2. 2-1)  can  be  looked  at,  in  general,  as  a  two  source  model  for 
the  vector  sequence  {u^},  whose  statistics  are  specified  by  the  statistics 
for  the  "source  inputs",  {e^  }  and{  b^}.  We  note  that  (3. 2. 2-1)  may  not 
always  qualify  to  be  a  vector  Markov  model 


by  defining 


Quj  -  Pu^,  -  fj 


f J  ’  EJ  +  "I  ‘ 


(3. 2. 2-3a) 


(3.2.2-3b) 
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-r—i 


This  is  because  the  sequence  {f^}  is  not  guaranteed  to  be  white,  even 
though  (e^)  i®  such.  From  the  linearity  of  (3. 2. 2-1),  the  two  source 
outputs  can  be  identified  by  the  solution  of  the  equations 


S  ■  S-1 +  s 


(3.2.2-4a) 

(3.2.2-4b) 


and  the  image  field  is  the  composition  [  31  ] 


o  b 

UJ  "  “J  j  * 


C3. 2.2-5) 


The  second  source  output  u^  depends  completely  on  the  boundary  variables. 
Hence,  as  the  sire  of  the  image  column  gets  larger,  one  would  expect  a 
smaller  contribution  of  u^  in  influencing  the  statistical  properties  of 
Uj.  The  N  x  M  random  field  will  be  stationary  when  the  statistics  of  the 
boundary  variables  are  consistent  with  those  of  {u  .  In  other  words, 

mJ 

the  case  of  stationary  statistics  is  only  a  special  case  of  (3. 2. 2-5). 

For  minimum  variance  representations,  the  orthogonality  condition 

(2. 3. 1-1)  implies  that  the  two  source  outputs  u°  and  u^  are  orthogonal. 

For  stationary,  white  noise  driven,  semlcausal  models,  one  does  not  have 

this  condition.  If  one  starts  by  specifying  the  boundary  variables  to 

be  orthogonal  to  the  stationary  white  noise  field  {e,  .) ,  then  the 

1 » J 

finite  N  *  M  field  {u  .}  need  not  be  stationary.  Asymptotically, 

mJ 

(N,M^“),  (3. 2. 2-4)  and  (3. 2. 2-5)  will  yield  stationary  random  field 

(u  }  (assuming  Q  is  always  nonsingular  and  (3. 2. 2-1)  is  a  stable 
1  *  J 

recursive  relation). 
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3.2.3  Stochastic  Decoupling  of  Semicausal  Models  by  the  Cosine  Transform 

for  a  large  number  of  image  fields,  the  adjacent  elements  of  a 
column  (or  a  row)  are  highly  correlated  and  often  the  correlation  param¬ 
eter  has  a  value  around  0.95.  Thus,  for  a  finite  N  *  M  random  field,  it 
is  reasonable  to  assume  that  the  boundary  variables  that  need  to  be 
specified  in  (3. 2. 2-1)  satisfy  the  conditions 


uo.j  “  ui,j 

Vl.j  “  “b.j  '  (3. 2. 3-1) 

This  means  the  two  outermost  boundary  rows  of  the  (N+2)  x  M  field  fu.  ; 

*•  •  J 

0  <_  i  <_  N+l,  1  <  J  <  M)  are  equal.  This  also  means  that  (3. 2. 2-1) 
subject  to  (3. 2. 3-1)  is  no  longer  a  stationary  field  on  the  (N+2)  *  M 
array  of  points.  Using  (3. 2. 3-1)  in  (3. 2. 3-1)  and  rearranging  terms, 
we  get  a  one  source  model  (the  second  source  is  deterministic,  as 
defined  by  (3.2. 3-1)  ) 


(3. 2. 3-2) 


(3. 2. 3-3.) 
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and 


r  i . 

for 

SCI 

phV 

for 

SC? 

(3.2.3-3b) 


For  SCI  model,  (r  }  Is  still  assumed  to  be  a  stationary  white  noise 
•  »  J 

field.  For  SC2  model,  (3.2. 3-2)  should  satisfy  the  minimum  variance 
conditions  of  (2. 3.1-1).  When  applied  to  (3. 2. 3-2),  this  Implies 


Kick  -  6Vj,v 


Hence,  the  covariance  of  is  now  given  as 


R  -  Ec.c‘ 
e  J  k 


B  I«.  .  ,  for  SCI 

J 


B  QcOjk,  for  SC2 


(3. 2. 3-4) 


Proposition 

The  KI.  transform  of  any  vector  ti^  of  the  random  field  of  (3. 2. 3-2) 
is  the  discrete  Cosine  transform  (DCT)  (  1  ]  for  SCI  as  well  as  SC2 
models.  For  proof,  see  (  35  ]. 

The  eigenvalues  of  are  given  by  (see  (4.2-llc)] 


X  .  ■  1  -  2o  cor 
cl 


1  <  i  <  N. 


(3. 2. 3-5) 


Define  ?c  to  be  the  DCT  matrix.  Since  it  diagonalizes  the  matrices  Qc 
and  P  we  have 

▼cQcY^  -  -  Diagonal  fA^)  (3.2.3-6a) 
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1 


T 

Y  P  w  * 

C  ■  c 


-  rc  ■  <V 


Y  ,  for  SCI 
p.  \  .  ,  for  SC2 

n  Cl 


(3.2.3-6b) 


Also,  define 


v,  •  Y  u  . 

j  c  J 


(3.2. 3-7) 


Multiplying  both  sides  of  (3. 2. 3-2)  by  ?  and  using  (3. 2. 3-6)  and 

c 

(3. 2. 3-7)  we  obtain 


where 


Acvj  •  rVi  +  'j 


et  ¥  c.  . 

J  c  j 


(3.2. 3- 8a) 


(3. 2. 3- 8b) 


Since  A  and  T  are  diagonal,  (3.2.3-8a)  is  a  set  of  decoupled  equations 


XCivJ<i)  *  Yivj-l(i)  +  ej(1)  1  <  1  <  N.  (3. 2. 3-9) 


From  (3.2. 3-4)  and  (3.2.3-6a)  and  (3.2.3-8b),  {c^(i)}  is  a  white  noise 
field  which  satisfies  the  equation 


EejCDe^U)  -  ■{ 


1  •torSC1 


(3.2.3-10) 


8  X  6  6.  ,,  for  SC2  . 

d  J  t  *  I*** 


Hence,  for  each  i,  the  sequence  (v^(i),  V  j)  is  a  first  order  Markov 
sequence.  Moreover,  if  the  initial  variables  {v0(i),  1  <  i  <  M)  are  un¬ 
correlated  then  for  V  j  >  0,  {v^(i),  1  <  i  <  N)  will  be  uncorrelated 
elements  of  the  vector  v  ^ .  Also,  for  •  0  or  u^  ■  u—  (i.e.,  steady 
state  in  J),  the  above  condition  will  always  be  satisfied.  For  arbitrary 
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initial  conditions,  say  -  c,  where  c  is  an  arbitrary  random  vector, 
Uj  will  have  a  decomposition 


'  UJ  *  “j 

(3.2.3-lla) 

where 

■  S-i +  ci 

c 

©  o 

1 

o 

(3.2.3-llb) 

and 

Qc“l 

■  S-1 

u0  -  c  . 

(3.2.3-llc) 

The  KL  transform  of  u°  now  will  be  the  DCT.  In  a  practical  DPCM  system, 
the  statistics  of  the  initial  condition  do  not  affect  the  ultimate 
system  performance  because  the  initial  vector  can  always  be  quantized 
accurately  and  transmitted  before  starting  the  DPCM  transmission. 


3.2.4  Stochastic  Decoupling  of  Stationary  Semicausal  Fields 

Random  fields  for  which  (3. 2. 3-1)  is  not  a  reasonable  assumption 
(i.e.,  one  step  correlation  in  the  "j"  direction  is  not  close  to  1)  or 
when  stationarity  of  the  random  field  is  an  essential  requirement  of  the 
model,  we  have  to  restrict  ourselves  to  (3. 2. 2-1).  The  boundary  variables 
(elements  of  b^)  are  random  variables  which  are  samples  from  the  station¬ 
ary  random  field  {u  )  and  cannot  be  specified  Independently.  In  that 
i » J 

case,  (3.2. 2-1)  admits  the  decomposition  (see  (3.2. 2-4)  -  (3. 2. 2-5)] 


S  •  S-i  *  cj 

e 

o  o 

■ 

o 

J 

>  1 

9u)  *  Pul-1  +  kJ 

b 

uo  ’  uo 

j 

>  1 

where 

o  .  b  .  .  , 

3  j  j  J  - 

and  uQ  is  the  Initial  value  of  the  sequence  (u.). 


(3.2.4-la) 

(3.2.4-lb) 

(3. 2.4-2) 
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3.3  Hybrid  Coder  Design  for  Seaicausal  Random  Fields 

For  the  sake  of  simplicity  of  presentation,  we  will  only  consider 
the  semlcausal  models  of  (3. 2. 3-2)  which  are  decoupled  by  the  Cosine 
transform.  The  sequence  of  decoupled  equations  (3. 2. 3-9)  and  (3.2.3-10) 
can  be  rewritten  as 


3.3.1  DPCM  Equations 

Since  (3.3-1)  is  a  sequence  of  N  independent  Markov  processes,  each 
of  these  sequences  could  be  transmitted  via  an  independent  DPCM  channel, 
as  shown  in  Fig.  3.2.  The  equation  at  the  transmitter  and  receiver  are 
simply  as  follows: 

Predictor:  v.j(i)  “  pivj-l^  (3.3.1-la) 

Differential  Signal:  e^d)  ■  v^(l)  -  Vj(i)  (3.3.1-lb) 
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(3.3.1-lc) 


O 


* 

Quantizer  Output:  e^d) 

Reconstruction  Filter;  v^(i)  »  p^v.  +  ej(D  (3.3.1-ld) 

The  encoding  schetne  requires,  first,  to  take  the  Cosine  transform  of  each 
column  vector  u^.  The  DCT  is  a  fast  transform  and  is  also  the  K.L  trans¬ 
form  of  each  column  of  the  semicausal  image  fields  characterized  by 
(3. 2. 3-2).  This  is  followed  by  N  DPCM  channels  for  predictive  coding  of 
successive  Cosine  transformed  vectors.  The  receiver  simply  reconstructs 
the  transformed  vectors  according  to  (3.3.1-ld)  and  performs  the  inverse 
Cosine  transformation.  To  complete  the  design,  we  now  need  to  specify 
the  quantizers  in  the  various  DPCM  channels. 


3.3.2  Bit  Allocation 

We  will  assume  all  the  quantizers  are  identical  in  their  charac¬ 
teristics  (i.e.,  bit  rate  vs.  distortion).  Let 


P  "  Average  desired  bit  rate  in  bits/pixel  of  the  image  field. 

f (n)  •  Mean  square  distortion  of  a  quantizer  with  2n  levels  for 

unit  variance  input  random  variables.  This  is  a  positive, 
monotonically  decreasing  convex  function.  We  will  assume 
f(n)  has  a  continuous  first  derivative. 

0  <  n^  ”  Number  of  bits  per  pixel,  allocated  to  the  ith  DPCM  channel. 

From  the  definition  of  P,  we  have 

1  N  ^ 

i  "l  ‘  P  "i  i  0  •  (3.3.2-U 


The  average  mean  square  distortion  in  encoding  of  the  image  vector  u  Is 
given  by 
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D  -  i  EiUjiiij 


(3. 3. 2-2) 


where 


Since  Vj  -  ^cuj  *s  a  un*tary  transformation  and  v^(i)  are  coded  via 


DPCM,  we  have 


6Vj(i)  -  Vj(i)  ~  Vj(i) 

■  vj<^  “  "  ej<*> 

-  e^i)  “  Pj(i) 


-  fle^O) 


(3. 3. 2-3) 


i.e.,  the  reconstruction  error  in  v^(i)  (assuming  noiseless  channel) 
equals  the  quantization  error  in  the  ith  channel.  This  gives 


D  ■  s  K''ISv) 


1  N  2 
-  jj  l  R(6e  (i))Z  . 
i-1  J 


(3. 3. 2-4) 


Thus,  the  average  distortion  per  pixel  is  the  average  quantizer  distor- 
2 

tion.  If  P^(i)  denotes  the  variance  of  the  jth  quantizer  input  in  the 
ith  DPCM  channel,  then 


EttejO))2  -  e2(i)f(ni)  ¥  J  . 


(3. 3. 2-5) 


Using  (3.3-1),  (3.3.1-la)  and  (3. 3.2-4)  in  (3.3.1-lb)  and  simplifying, 
we  get 
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(3. 3. 2-6) 


e^d)  -  Vj(i)  -  v.  (i) 

-  (i)  +  Pj ' e  j_j  • 

Eqns .  (3. 3. 2-3)  and  (3. 3. 2-6)  together  give  a  recursive  relation  for  the 
variance  of  c^O)  rts 

Bj(l)  -  B2(l)  +  §0  -  0  (3. 3. 2-7) 

where  we  have  used  the  fact  that  c^(i)  is  uncorrelated  with  the  quanti¬ 
zation  error  at  step  J-l.  In  steady  state 

B2(i)  ‘  0 2 < i ) - -  .  (3.3. 2-8) 

j  1  -  pjftn^ 

Hence,  (3. 3. 2-3)  becomes 

.  (i)f(n.) 

E(&e  (i))  "  5  (3. 3. 2-9) 

j  1  -  Pjf(ni) 


giving 


D  - 


1 

N 


(3.3.2-10a) 


where 

f(x) 

g  (x)  "  - « -  x  >  0  .  (3.3.2-10b) 

1  -  Pjftx) 


The  function  gj(x)  represents  the  mean  square  distortion  of  the  ith  DPOi 
loop  of  a  Markov  process  whose  prediction  error  variance  is  unity.  Like 
the  function  f(x),  for  |pj  <  1,  each  (x)  is  also  a  positive,  monotoni- 
cally  decreasing,  convex  function  with  a  continuous  first  derivative. 
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The  bit  allocation  problem  la  to  minimize  (3.3.2-10a)  subject  to 
the  constraints  of  (3. 3.2-1).  This  is  equivalent  to  finding 

1  N  ■> 

D'  -  min -.XP  -  -  \  [g,(n  )fl  (i)  +  Xn  ]*) 
nt>0  i-1  1  1  11 

X 


N 

•  min\Xp  -  —  \  max  (gi(ni)SZ(i)  +  Xn^  •  (3.3.2-11) 


I 


i-1  ni>o 


Solution  of  constrained  minimization  problems  of  this  type  are  well  known 
in  optimization  theory.  Recently,  Segall  [  66  j  has  considered  a  similar 
bit  allocat ion  problem  [wi th  g^(x)  -  g(x),  V  i]  for  encoding  of  vector 
sources.  The  solution  is  given  by 


gi.(ni>®2(!)  +  X  "  0  A  <  -g;<0)B2(i) 


"i  "  ° 


X  >  -g’(O)0*(i)  . 


(3.3.2-12) 


Deno  t ing 


h^x)  “  g[  1(*) 


(3.3.2-13) 


we  can  solve  (3.3.2-12)  for  as 


.  X  <  -gj(0)8  (1) 


"i  "  ' 

0  .  X  >  -g!(0)B2(i) 

1 

where  X  is  the  root  of  the  nonlinear  equation 


(3.3.2-14) 


l  M- 

i:X<-gj(o)e2(i) 


B2(i) 


NP 


(3.3.2-15) 
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Th«?  value  of  the  minimum  distortion  is 


where  0  Is  the  solution  of 


I  max  [0,  Ij1or2  - 


(3. 3 . 2— 21b) 


and  the  minimized  distortion  Is 


“..n-S  wW) 


(3.3.2-21c) 


Eqns.  (3.3.2-21)  arc  valid  for  both  SCI  and  SC2  models,  where  B  (1)  is 
defined  by  (3.3-2b). 


2.  Consider  the  case  when  f(*)  and  g^(*)  are  given  by  (3.3.2-18)  and 


(3.3.2-10b)  respectively.  Then 


*;<*> 


f'(x) 

(1  -  p2f(x))2 


(-21ogc2)2  A 
(1  -  p22"2x)2 


(3.3.2-22) 


From  this  equation  one  can  obtain  a  unique  inverse  (note  that  g^(x) 


must  be  nonpositive  for  all  x  >  0) ,  ns  (see  Appendix  C) 


-Slog 


t!  1 (x)  » h, (x) 


E-1or  2+P2x+/(1or  2)1 

_ e _ i  « 

4 

P.x 


-  (21or^2)o* x 


-S1or2 


21orc2 


,  Pj  -  0  . 


■  *  °t4 


(3.3.2-23) 


For  SCI  model,  defining 


(3.3.2-24) 
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and  using  (3.3.2-23)  In  (3.3.2-14),  (3.3.2-15)  and  (3.3.2-16)  we 


t lm 


Slog,  w(i ,0)  ,  0  <  0  <  -  ^!r-= 

‘0  (1  -  P  y 


,  e  >  - 


2,2 

(1  -  Pt) 


(3.3.2-25a) 


where 


w(i,0)  $ 


1  .  Pt  -  0 


i  r. .  2pi°  i.  Api° 

,4  1  +  «2~  ‘  l+02/jV  *  P1 


2pi  L  8  (0  \  8  (l) 


3.2-25b) 


fl  Is  the  solution  of 


l  Slog2  -  NP 

i  :0<6<  -£-4^5 

(1-Pj) 


(3.3.2-25c) 


and  the  minimum  distortion  Is 


Dmin  "  h'  £  «in(82(i),  G/w(i,9)]  . 


(3. 3. 2-25d) 


In  the  case  of  SC2  model,  substitution  of  p  ■  p  •  constant, 

1  n 

(3.3.2-25)  above  gives  the  desired  bit  allocation  and  distortion 
formulas.  The  bit  allocations  are  obtained  by  first  solving 
(3.3.2-25r)  iteratively  (by  Newton's  method,  for  example)  to  obtain 
0,  for  a  fixed  rate  P.  Given  0,  (3.3.2-25a)  and  (3.3.2-25b)  yield 
the  desired  allocations  {n^}  and  (3.3.2-25d)  gives  the  corresponding 
distortion. 
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3.3.3  Integer  Bit  Allocation 

In  the  above,  we  assumed  to  be  a  real  number.  In  many  practical 
situations  n^  is  required  to  be  an  Integer.  An  approximate  integer  bit 
allocation  can  be  obtained  via  the  foregoing  method  by  rounding  n^  to  the 
nearest  integer.  The  actual  bit  rate  then  would  be 


\  "s  l  <"i> 

i-1 


where  [x]  indicates  the  nearest  integer  to  x.  To  obtain  the  exact  solu¬ 
tion,  the  minimization  of  (3.3.2-lOa)  is  subject  to  the  constraint  that 
{n^}  be  a  set  of  nonnegative  integers.  The  solution  is  obtained  via  an 
iterative  algorithm  due  to  Fox  [  11  ]  as  follows. 


Let  n  denote  a 


vector  of  elements  fn^,  1  <  i  <  ll}  . 


1.  Start  with  the  allocation  n  -  0. 


Set  1*1. 


3.  n*'  *■  n^  +  e^  where  e^  is  the  ith  unit  vector  and  is  any  index  for 
which 

Vi "  e2(i)*i(n[’1)  -  +  i) 

is  maximum.  The  quantity  A  .  is  called  the  marginal  return  of 
increasing  the  1th  allocation  in  Ith  step  by  one  bit. 

H  t 

4.  If  £  n*-  >  PN,  terminate;  otherwise  (  *1+  1  and  go  to  step  3. 

i-1 

This  algorithm  simply  means  the  marginal  returns  A.  .  are  calculated 

*  t  *■ 

for  1  <  1  <  »  and  i  *  1,2,...,  and  are  arranged  in  a  decreasing  order 
until  all  the  bits  arc  exhausted.  Tables  3. land  3.2  show  the  bit 


allocations  for  the  SCI,  SC2  ami  separable  covariance  models  for  Girl 
image  obtained  via  real  bit  allocation  and  integer  bit  allocation  algo¬ 
rithms  described  above. 

3.4  Adaptive  Hybrid  Coding 

The  coding  scheme  of  the  previous  section  can  be  adapted  to  image 
fields  whose  spatial  statistics  vary  slowly  by  updating  the  various 
coefficients  of  the  image  model  with  variations  of  the  statistics.  A 
complete  update  of  the  image  model  parameters  could  often  increase  the 
complexity  of  the  coder  to  make  it  impractical.  In  the  sequel  we  con¬ 
sider  two  adaptive  schemes  which  offer  reasonable  improvement  in  coder 
performance  with  only  a  marginal  increase  in  its  complexity. 

3.4.1  Adaptive  Variance  Estimation 

For  a  fixed  predictor  in  the  feedback  loop  of  a  DPCM  channel,  the 
variance  of  the  prediction  error  will  fluctuate  with  changes  in  image 
statistics.  In  this  scheme  we  update  the  variance  of  the  prediction 
error  at  each  step  j.  This  updated  variance  is  used  to  adjust  the 
sparing  of  the  quantiser  levels  in  each  DPCM  channel.  For  mean  square 
error  criterion,  the  decision  and  reconstruction  levels  of  a  quantizer 
are  directly  proportional  to  the  standard  deviation  of  the  Input  random 
variables.  Hence,  the  above  mentioned  adaptation  can  be  achieved  by 
aimply  normalizing  the  prediction  error  by  its  updated  standard  deviation. 
The  quantizer  levels  are  then  fixed  and  correspond  to  unit  variance  ran¬ 
dom  variables.  Fig.  3.4  shows  this  scheme. 
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Let 


oj(i)  ^  Variance  of  e^(i),  the  prediction  error  at  step  j  of  the 


ith  DPCM  loop  -  E{e  (i)]' 


(3.4.1-la) 


ItO  A  t  0 

o^(i)  “  Variance  of  the  quantised  values  “  E [ e^ (i ) )  -  (3.4.1-lb) 


Since  the  quantised  variables  (i)  are  available  both  at  the  receiver  and 

*2 

the  transmitter,  it  is  easy  to  estimate  0^(1).  Caatellino  et  al  [  *  ] 

*2 

have  suggested  a  formula  for  an  estimate  of  0^(i)  (for  single  channel 
DPCM  system)  as 

J 

$?(i)-^*  l  Y  e*2..(i)  +  Y3SJ(i)  0  <  y  <  1.  (3.4. 1-2) 

J  Y  »  J  -mri  1 


This  has  been  called  an  exponential  average  variance  estimator.  A  more 
convenient  form  of  (3.4. 1-2)  is  the  recursion 


^+1(i)  -  (l-y)e*2(i)  +  Y*j(i)  j  -  1,2,... 


(3.4. 1-3) 


For  small  quantization  errors  one  may  take  to  he  equal  to  o^.  For 
Lloyd-Max  quantizers  [  44  1  a  more  accurate  estimate  of  is  possible. 
Since  the  variance  of  a  Lloyd-Max  quantizer  input  equals  the  sum  of  the 
variances  of  the  quantizer  output  and  the  quantization  error,  we  have 


A  2 

•j(» 


A  2 

0j(i) 


Sj(i)  +  f(nt)Sj(i) 


fro 

l-f(n4) 


ni  >  1 


(3. 4. 1-4) 


where  is  the  number  of  bits  allocated  to  the  ith  channel.  Using 
(3.4. 1-4)  in  (3.4. 1-3)  we  get 
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°j+l(i>  "  idfOO  ej2(i)  +  ySj(1)  ni  >  1  0  <  Y  <  1.  (3. 4. 1-5) 


The  above  estimate  become  poor  for  DPCM  channels  which  are  assigned  small 

/V 

number  of  bits  (n^  ~  1).  For  these  channels  0^(1)  could  be  estimated  by 

A 

some  sort  of  extrapolation  of  the  already  estimated  0^(1).  For  example, 
for  the  SCI  model ,  we  know 


0^(1)  *  E  %2Al) - j  - - 5~ 

3  j  [l-pjf(ni)]x2i 


i  e  7 


(3.4. 1-6) 


where  7  Is  the  set  of,  say,  m  channels  for  which  n^  >  1.  From  (3. 4. 1-6) 
2  -  2 

an  estimate  of  6  at  step  ),  denoted  by  P^  is  given  by 


“j  ■  S  j,3I<i,xo'1  -  ^f<V  • 


(3.4. 1-7) 


This  gives  the  extrapolated  estimates  as 


-2  J 

0«(D  -  ~2 - 2 - 

i  Xc1Il‘°if(n1)1 


i  i  7 


(3. 4. 1-8) 


For  SC2  model,  the  equations  corresponding  to  (3.4. 1-7)  and  (3.4. 1-8) 
become 


’  ;  j, 


(3.4. 1-9) 


A  2 

Oj(i) 


xci*Kf(ni)3 


i  i  l  . 


(3.4.1-10) 


The  adaptive  variance  estimator  only  needs  to  solve  Eqns.  (3. 4. 1-5), 
(3. 4. 1-7)  and  (3. 4. 1-8)  for  SCI  model  and  (3. 4. 1-5),  (3.4. 1-9)  and 


0 
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(3.4.1-10)  for  SC2  model.  Ue  finally  note  that  this  adaptive  scheme 
maintains  the  same  constant  data  rate  in  each  DPCM  channel  as  in  the 
nonadaptlve  scheme  of  the  previous  section. 


3.4.2  Adaptive  Classif tear  ion 


In  this  method  each  image  column  is  classified  as  belonging  to  one 
of  K  predetermined  classes.  Each  of  these  classes  is  determined  accord¬ 
ing  to  the  activity  in  that  image  column.  The  activity  in  an  image  column 
is  measured  by  the  variance  of  that  column.  For  a  given  class  of  Images 
the  probability  distribution  function  of  the  variances  of  all  the  columns 
can  be  pre-determined.  (Note  that  this  probability  distribution  function 
could  depend  on  the  size  of  the  image  column  used  in  hybrid  coding.) 

In  the  DPCM  loop,  the  transformed  signal  in  each  DPCM  channel  is 
processed  as  before  except  that  the  differential  signal  is  quantized 
according  to  the  classification  of  that  image  column.  Thus,  the  quantizer 
in  the  DPCM  loop  is  adapted  to  the  image  activity  in  each  column.  The 
number  of  quantization  levels  (or  bits)  depend  on  the  variance  of  a  class 
so  that  columns  of  high  dynamic  activity  are  assigned  more  bits  than 
those  of  low  activity.  In  this  way  quantizing  bits  are  employed  more 
efficiently.  However,  the  classification  information  for  each  column 


needs  to  be  transmitted.  This  can  be  done  by  sending  an  extra  logjK  bits 
per  column  of  ~  logJC  bits/pixel.  Fig.  3.5  shows  the  adaptive  classi- 
f leaf  ion  scheme.  Experimentally  it  was  found  that  the  variance  of  suc¬ 
cessive  columns,  denoted  by  {£j}  (say),  are  highly  correlated.  Instead 
of  transmitting  the  classif ication  map,  a  DPCM  loop  for  transmission  of 


are  used  for  classification  at  the  receiver  as  well  as  the  trans¬ 


mitter.  A  two  bit  quantizer  was  found  to  give  very  accurate  reproduction* 


Let  p(0  ■  probability  density  of  the  variances  £  of  a  column. 


Suppose  the  K  classes  have  been  predetermined  and  let  (t^.t^^). 

1  <_  Ic  <_  K  denote  the  decision  boundaries  for  the  kth  class.  Since 
C  >  0,  -  0,  e... .  -  Then 


p(£)d£  -  probability  of  the  kth  class 


(3.4. 2-1) 


1  f  k  1 

T  "  —  £p(£)d r,  -  variance  of  the  kth  class. 

k  Pk  ) 

*k 


(3.4. 2-2) 


First  assume  that  the  thresholds  { t^)  has  been  somehow  predetermined  and 

p,  and  C.  are  known.  Also,  we  will  assume  that  for  each  class  k,  the 
k  k 

image  model  parameters  have  been  predetermined  and  are  known.  To  this 

2 

end,  let  us  denote  by  p  ,  8  (i.k)  n  ,  etc.,  the  quantities  for  the 

1  »  K  1  |  K 

ith  DPCM  channel  and  kth  class  corresponding  to  the  definitions  of  p^, 

2 

8  (i),  n^  etc.,  respectively.  Then,  following  the  development  of  the 
foregoing  sections,  we  can  write  the  expected  average  distortion  for  any 
image  column  (in  steady  state  DPCM)  as 


w  K 

D  “  N  J  i  B  <i,k)gi,k(ni,k)pk 

1 "1  k“ 1 


(3.4. 2-3) 


where  n  .  >  0  are  such  that  the  expected  average  bit  rate  must  be  con- 
1 1  k 


slant,  ie.. 


N  K 


N  l  l  ni,kPk  "  ?  * 

i-1  k»l 


(3.4. 2-4) 


The  bit  allocation  problem  associated  with  (3.4. 2-3)  and  (3. 4. 2-4)  can 
be  solved  along  the  same  lines  as  in  section  3.3.  As  an  example,  consider 
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the  case  when 


8i  *  f<x>  "  2_2X 

then  one  obtains 

"l.k  "  Bax  |^0»  **lQg2~  0.4. 2-5a) 

where  0  is  the  solution  of 

l  l  wax  [ 0,  Slog/-  y*  p.  -HP  (3.4.2-5b) 

i-1  k-1  L  0  J  k 

and  the  mlnLmua  achievable  distortion  is 

N  K 

D  ,  -  l  l  «nln[0,  62(i.k)]p.  .  (3.4.2-5c) 

“■  1-1  k-1  k 

This  method  of  adaptive  hybrid  coding  requires  1)  measurement  of  variance 
of  each  column  and  its  transmission  by  an  extra  DPCM  channel,  li)  classi¬ 
fication  of  each  column  to  one  of  K  classes  based  on  thresholds  {t^}, 
and  ill)  storing  of  K  bit  allocation  tables,  one  for  each  class. 

3.5  Experimental  Results 

Several  computer  experiments  have  been  performed  on  the  Girl  and 
Chemical  Plant  images  to  simulate  the  hybrid  coding  schemes  discussed  in 
the  previous  sections.  Table  3.3  shows  the  parameters  of  the  various  image 
models  used  here.  In  all  the  experiments,  the  Cosine  transform  based 
semicausal  models  were  used.  A  256  x  256  image  was  divided  into  16  strips, 
each  of  sise  16  x  256.  Each  strip  was  hybrid  coded  Independently.  This 
way  only  a  16  step  Cosine  transform  Is  needed  which  makes  the  scheme 
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feasible  for  a  real  time  implementation  [  46  1  •  In  all  the  performance 
measurements,  the  signal  to  noise  ratio  (SNR)  is  defined  as 


SNR  -  201or10  ( 


Peak  to  Peak  Value  of  the  Signal^. 


r.m.s.  error 


(3.5-1) 


For  an  eight  bit  original  data  this  becomes 


SNR  -  201 og 


10  r.m.s.  error 


(3.5-2) 


The  normalized  mean  square  error  is  measured  over  the  entire  256  x  256 


image  and  is  defined  as 


256  256 


*  \  2 


l  l  (ui  j  ~  ui 
K>,SF  -  i-JU-i—  - _ _ 

l  l  U1  4 

i-1  j-1  l,J 


(3.5-3) 


In  each  DPCM  channel  a  compandor  is  used  as  an  approximation  to  the 
Max  quantizer  (see  Appendix  A). 


3.5.1  Nonadaptlve  Hybrid  Coding 

Figure  3.6  shows  the  SNR  vs.  bit  rate  and  rate  distortion  curve  for 
the  SCI,  SC2  and  the  separable  covariance  image  models  for  the  girl 
image.  As  expected,  the  performances  of  SC2  and  separable  covariance 
models  arc  very  close.  The  SCI  model  performs  better  than  the  other  two 
by  about  1.5  to  2  dB,  in  general,  above  rates  of  1  bit/pixel.  The  same 
set  of  plots  for  Chemical  Plant  image  is  shown  in  Fig.  3.7.  Since  the 
performances  of  SC2  and  separable  covariance  models  are  close,  for 
simplicity,  the  curve  for  the  separable  covariance  model  is  not  shown  in 
Fig.  3.7. 


3.5.2  Adaptive  Variance  Estimation 


Figure  3.8  shows  the  SNR  vs.  bit  rate  and  rate  distortion  curve 
for  the  various  models  of  the  Girl  Iamge  for  this  case.  In  these 
experiments,  the  variance  estimates  were  updated  only  for  channels  for 
which  n^  2  bits/pixel. 

The  performance  of  SC2  and  separable  covariance  model  improves  by 
about  1.5  dB  at  1  bit  rate  and  about  4  dB  at  2  bit  rate.  These  values 
can  be  found  from  Table  3.4,  which  gives  the  quantitative  comparisons 
on  a  model  by  model  basis.  The  reason  for  larger  improvment  at  2  bit 
rate  is  that  the  variance  estimates  (based  on  quantised  data)  arc  more 
accurate  since  many  more  channels  have  n^  >  2  bits/pixel.  The  per¬ 
formance  improvement  for  SCI  model  is  less  (only  0.5  dB  at  1  bit  rate 
and  1.5  dB  at  2  bit  rate).  This  could  he  because  the  SCI  model  per¬ 
formance  was  already  better  (than  SC2)  to  start  with  and  the  margin  of 
improvement  by  updating  prediction  error  variance  is  small.  As  the  size 
of  the  image  column  (N)  is  increased,  the  performance  of  the  adaptive 
variance  estimator  was  found  to  improve.  This  is,  again,  because  the 
number  of  channels  (for  a  fixed  average  rate)  with  bit  allocation  n^  >.  2 
would  increase  giving  a  better  estimate  of  the  variances  (as  well  as  of 
the  extrapolated  values).  Similar  conclusions  can  be  made  for  the 
Chemical  Plant  Image.  The  corresponding  results  are  shown  in  Fig.  3.9 
and  Table  3.5. 

3.5.3  Adaptive  Classification 

Two  different  experiments  utilizing  this  method  were  performed.  Tn 
the  first  experiment,  a  probability  distribution  function  of  the  column 
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variances  was  determined  from  the  histogram  of  their  sample  variances. 


Classification  thresholds  { t^)  were  determined  for  a  K  ■>  4  case  assuming 
equal  probability  for  each  class.  Then  each  column  of  the  image  was 
classified.  The  actual  measurements  model  for  the  Cosine  transformed 
vectors  was  assumed,  as 


v^(i,k)  -  (i,k)  +  (i  »h) 


1  <  i  <  N  .  k  -  1....4  (3. 5. 3-1) 


wi  th 


Eej ( i ,k)  B2(i,k) 

-  (l-P^k)BfM  (i.k)]2  . 


(3. 5. 3-2) 


The  parameters  p^  ^  and  B  (i,k)  were  measured  for  each  class  and  bit 
allocation  wa3  done  according  to  the  method  described  in  the  previous 
chapter.  Fig.  3.5  shows  the  adaptive  classification  hybrid  coding 
scheme.  Fig.  3.10  shows  the  16  *  32  classification  map  of  16  x  256 
portion  of  the  image  with  16  element  columns.  (Recall  that  the  given 
256  x  256  imago  was  divided  into  sixteen  16  *  256,  rectangular  subimages.) 

This  method  gave  a  large  improvement  (see  Fig.  3.11  for  the  Girl 
Inage  and  Fig.  3.12  for  the  Chemical  Plant  Image)  in  SNR  over  the  pre¬ 
vious  two  schemes  (about  5.5  dB  at  1  bit  rate  and  8  dB  at  2  bit  rate 
over  nonadapt ive) .  However,  the  above  procedure  is  somewhat  impractical 
for  online  transmission  (could  be  acceptable  for  storage)  because  it 

A 

requires  identification  of  the  2NK  -  8n  parameters,  p  .  and  E(v,  (l,k)  *"  ] 

J 

for  each  image.  These  parameters  are  then  used  to  determine  the  bit 
« I  locations  for  various  classlf ications.  Assuming  eight  bits  are  re- 
■  '  to  transmit  each  of  these  parameters  (this  is  actually  a  pessi- 
.  »s  1 1  mate) ,  the  additional  overhead  for  their  transmission  is  only 
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2  2 
where  is  the  variance  of  the  columns  belonging  to  the  kth  class,  8^(1) 

represents  the  prediction  error  variance  of  the  ith,  unit  variance  Markov 
process.  Compared  to  the  nonadaptive  case,  the  only  additional  computa¬ 
tion  this  model  requires  is  the  measurement  of  the  sample  variance  of  an 
image  column  and  classifying  it.  Figs.  3.13  and  3.14  show  the  bit  allo¬ 
cations  for  various  models  used  for  adaptive  classification  hybrid  coding. 
Since  the  first  DPCM  channel  in  the  Cosine  transform  represents  the  d.c. 
value  (or  sample  wean  of  the  columns),  the  number  of  bits  allocated  to 
this  channel  (for  a  given  bit  rate  and  class)  were  kept  the  same  for  all 
different  models.  The  bit  allocation  algorithm  of  previous  sections  was 
thus  used  only  for  i  ■  2,...,Nlnthe  case  of  SCI  and  SC2  models.  Figs. 
3.11  and  3.12  show  the  performance  of  these  models  when  used  for  adaptive 
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classif ication.  Generally,  the  performance  of  the  SCI  model  Is 
superior  to  that  of  the  S C2.  For  the  Girl  Image,  It  Is  better  than 
nonadaptive  (separable  covariance  model)  hybrid  coding  by  about  3  dB, 
but  is  worse  than  the  previous  adaptive  classification  experiment  by 
2  to  3  dB.  This,  of  course,  is  the  price  the  designer  has  to  pay  for  a 
simplified  system.  However,  for  the  Chemical  Plant  image,  the  adaptive 
class  if icatlon  scheme  does  not  give  much  SNR  improvement  over  the  non¬ 
adaptive  scheme,  it  is  worse  than  the  adaptive  variance  scheme.  This 
is  understood  by  referring  to  Figs.  3.15  and  3.16,  which  are  the  histo¬ 
grams  of  image  column  variances.  For  the  Girl  image,  the  variances 
have  a  large  dynamical  range,  whereas  the  variance  of  the  Chemical 
Plant  image  are  nearly  equal.  Hence  image  column  activities  of  the 
Chemical  Plant  image  remain  invariant  and  therefore  classification 
does  not  help.  Finally,  note  that  this  scheme  will  produce  variable  bit 
rate  from  one  image  column  to  the  next  because  of  a  possible  change  in 
classlf icat Ion. 

Figs.  3.17-3.19  show  the  results  of  various  hybrid  coding  schemes 
on  a  node]  by  model  basis  for  the  Girl  image.  Figs.  3.20  and  3.21  are 
those  for  the  Chemical  Plant  image.  Tables  3.4  and  3.5  give  the  quanti¬ 
tative  comparisons  of  the  above  mentioned  results.  Figs.  3.22  and  3.23 
show  the  encoded  Girl  image  at  approximate  bit  rates  of  1  and  2  bits/ 
pixel,  and  the  corresponding  figures  for  the  encoded  Chemical  Plant 
image  arc  shown  in  Figs.  3.24  and  3.25.  Fig.  3.26  compares  the  original 
and  the  adaptive  classified  (based  on  actual  measurement  model)  hybrid 
encoded  images  at  1  and  2  bits/pixel. 


Finally,  observe  that  as  an  image  column  size  is  Increased,  the 
ambiguity  in  classifying  it  into  different  classes  will  also  go  up 
because  the  activity  in  a  large  image  column  is  not  likely  to  be  uniform 
so  as  to  be  represented  by  its  variance. 

3.6  Extensions  of  Hybrid  Coding 

Hybrid  coding  techniques  can  be  generalized  and  extended  to  higher 
order  semicau sal  models  and  for  coding  of  noisy  images.  In  this  section 
ve  briefly  discuss  the  formulations  and  simulation  results. 

3.6.1  Hybrid  Encoded  Image  with  Channel  Error 

Since  practical  communication  channels  are  noisy,  it  ia  worthwhile 
to  study  the  channel  noise  effects  on  the  Hybrid  encoded  images.  A 
binary  syimnetric  channel  is  used  here  as  the  noisy  communication  channel. 
Fig.  3.27  shows  the  representation  of  such  a  channel,  where  the  probability 
of  receiving  an  Incorrect  symbol  is  and  a  correct  one  is  1  -  P^, 
regardless  of  which  symbol  is  transmitted. 

Fig.  3.28  shows  the  effect  of  channel  errors  on  hybrid  coding 
schemes.  Generally,  the  adaptive  classification  is  quite  robust.  The 
nonadaptlve  and  adaptive  variance  schemes  can  be  made  leas  sensitive  to 
channel  errors  by  reducing  somewhat  the  horizontal  correlation  parameter. 
The  SCI  model  adaptive  variance  estimation  aeemm  to  be  leas  sensitive 
than  the  SC2  model.  F'  s.  3.29  and  3.30  ahow  the  encoded  images  in  the 
presence  of  channel  errors.  Table  3.6  gives  the  quantitative  compari¬ 


sons  of  these  results. 


3.6.2  Higher  Order  Predictor  in  DPCM  Loop  of  Hybrid  Coding 

This  experiment  is  designed  to  have  a  nth  order  predictor  instead 
of  a  first  order  predictor  in  the  DPCM  loop.  Tabic  3.7  shows  the  results 
of  using  various  predictors.  Table  3.8  shows  the  numerical  values  of  the 
predictor  coefficients  which  are  precalculated  on  the  basis  of  the  statis¬ 
tics  of  the  transformed  image  strips. 

These  coefficients,  n^,  appear  in  the  nth  order  autoregressive 
process  (one  for  each  channel' 


•i  ‘  ^Vj-i  *  "j 


(3. 6. 2-1) 


and  are  evaluated  by  solving  the  linear,  Toeplitz  equations 


rk '  J/ k-i°i 


(3.6. 2-2) 


where  r^  ■■  Ev^v,^  is  the  covariance  of  the  ith  channel  process  in 


Hybrid  coding. 


The  variance  of  the  differential  signal  e^  is 


2  2  r 

°i  "  0  ■  i  Vi 

J  i-1  1  1 


(3. 6. 2-3) 


where  o  is  the  variance  of  the  input  image. 

A  block  diagram  of  this  scheme  is  shown  in  Fig.  3.31  and  Table  3.7 
shows  clearly  that  first  order  predictor  is  a  reasonable  choice  for 


Hybrid  Coding,  since  the  signal  to  noise  ratio  improvement  is  small  even 
though  the  order  of  the  predictor  goes  up  to  four. 
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3.6.3  Hybrid  Coding  of  Noisy  Images 

In  many  inage  transmission  applications,  the  source  image  may 
happen  to  be  noisy.  This  is  because  of  the  image  spot  pick-up  sensor 
characteristics  or  due  to  some  Inherent  limitations  inside  the  mechanism. 
The  foregoing  developed  coding  methods  can  be  modified  to  encode  noisy 
images,  if  ve  know  the  statistics  of  noisy  images  beforehand.  Let  us 
assume  a  given  noisy  image  degraded  by  additive  white  Gaussian  noise 
with  known  signal  to  noise  ratio  (S/N)  a  priori,  where  S/N  ia  defined  as 


S/N  - 


S tandard  deviation  of  the  noise-free  signal 
Standard  deviation  of  the  adding  noise 


The  additive  noise  in  the  given  image  can  be  filtered  by  replacing  a  pre¬ 
dictor  in  the  DPCM  loop  of  Hybrid  coder  by  a  discrete  Kalman  filter. 

This  design  allows  one  to  do  the  Joint  filtering  and  coding  simultan¬ 
eously.  Fig.  3.32  is  the  block  diagram  of  this  scheme.  The  message  and 

observation  models  in  the  ith  channel  of  Hybrid  coder  after  the  trans¬ 

formation  are  as  follows. 


Message : 


V j ( 1)  -  Ph(i)v  x(i)  +  tj(i) 


(3. 6. 3-1) 


Observation: 


*j(l)  -  v^i)  +  nj*^ 


(3.6. 3-2) 


-  °J<l-»h«»Vl,0Sk,0 


(3.6. 3-3) 


Eynyi+k)  -  (s/9)\06k0 


(3.6. 3-4) 


|m-n| 


{Aj)  are  the  diagonal  entries  of  eigenvalue  of  Ry  •  (pj^  ),  P^  and  p 
are  the  horizontal  and  vertical  corrections  respectively. 
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The  encoding  scheme  requires  first  taking  the  Cosine  transform  of  each 
noisy  Image  column  u^ ,  then  applying  a  bank  of  discrete  Kalman  filter 
DPCM  on  each  row  of  the  resulting  image.  Other  necessary  equations  used 
in  this  scheme  are  as  follows. 


Predicted  Observation: 

*j/j-l(1)  *  Ph(1)Vj-l(i) 

(3.6. 3-5) 

Innovation: 

( i)  “  Zj(i)  -  Zj /j— 1 ( i) 

(3.6. 3-6) 

Quantiser  Output: 

ej(i) 

(3.6. 3-7) 

Current  Estimate: 

VjO)  -  Ph(i)v  (i)  +  k.j(e)eJ(i) 

(3.6. 3-8) 

The  above  equations  are  very  similar  to  the  DPCM  equations  (page  46). 
Furthermore,  k^(i)  appears  in  (3. 6. 3-8),  is  called  the  Kalman  gain  and 
often  referred  to  as  a  correction  factor  to  the  innovation  e^d),  (it  is 
e^d)  essentially  in  this  scheme)  and  can  be  obtained  recursively  as 
follows. 


Cain : 


k  .  ( i)  ■  V-  (OlV-  (i)+V  (i))' 

J  J/J-1  J/J-l  J 


(3. 6. 3-9) 


A  priori  Variance : 


V-  (i)  -  p  (i)V-  (i)p. (i)  +  ▼  (3.6.3-10) 

Vl/J  h  Vj  h  € 


A  posteriori  Variance:  V-  (i)  ■  [l-k.(i))V-  (i) 

VJ  J  Vj/J-1 


(3.6.3-11) 


where 


Vl/J (1)  "Vi(i)  -  Vl/J (1) 


Vj(i)  -  V^i)  -  v  (i) 


(3.6.3-12) 

(3.6.3-13) 
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(3.6.3-14) 
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-  10 log  +  201ogl()  |  +  SNR1 

0 
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where  Is  the  variance  of  noise-free  Girl  Image  (In  our  example  it 
equals  1816.21).  Therefore,  the  above  equation  becomes 


SNR2  -  lOlog  ffffLj  ♦  201og10  |  *  SNR, 


(3.6.3-16) 


In  [34] it  has  shown  that  semicausal  image  models  can  be  used  to 
formulate  and  solve  problems  in  Image  restoration.  Here  various  hybrid 
coding  schemes  based  on  these  models  have  been  extended  to  encode  noisy 
images.  Due  to  the  knowledge  of  image  models,  hybrid  coding  schemes  can 
be  successfully  applied  to  encode  both  noisy  and  noiseless  images.  The 
results  of  this  section  promise  the  effectiveness  of  previously  developed 
hybrid  coding  schemes  in  image  data  compression  applications. 


3.7  Summary  of  Results  and  Conclusions 

Based  on  the  above  experiments  and  the  theory  discussed  in  this 
chapter,  the  following  conclusions  are  made. 

i)  It  was  found  that  the  performance  of  conventional  nonadaptive 
hybrid  coding  can  he  improved  up  to  5.5  dB  at  low  bit  rates 
(*1  bit/pixel)  by  improved  statistical  characterization  of  the  image 
via  more  accurate  models  and/or  via  adapting  the  parameters  of  the 
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chosen  nodel  to  changes  In  image  statistics.  The  actual  measure¬ 
ments  model  applied  in  the  adaptive  classification  coding  scheme 
gave  the  best  results.  However,  this  model  requires  identification 
of  a  large  number  of  parameters.  In  general,  the  SCI  model  per¬ 
formed  better  than  SC2  and  separable  covariance  models  by  2  to  3  dB. 

As  expected  from  theoretical  considerations,  the  SC2  and  separable 
covariance  models  performed  almost  identically. 

11)  The  adaptive  variance  estimation  technique  is  useful  when  it  is 

desired  to  maintain  the  performance  of  a  hybrid  coder,  designed  for 
a  nominal  statistics,  in  the  face  of  changing  statistics  (e.g., 
the  correlation  parameter).  If  the  Image  data  statistics  is  close 
to  the  design  parameters,  then  the  improvement  over  nonadaptive 
scheme  is  marginal  (1  to  1.5  dB  at  1  bit /pixel).  This  technique  is 
effective  at  relatively  higher  bit  rates  1  blt/plxel  and  for 
relatively  larger  image  column  sizes.  It  should  be  considered  not 
so  much  for  achieving  low  transmission  rates  as  for  achieving  higher 
SNR  at  moderate  rates.  The  complexity  over  a  nonadaptive  scheme  Is 
marginal,  and  It  can  be  hardware  Implemented  very  easily.  Finally, 
it  maintains  a  constant  bit  rate  from  one  image  column  to  the  next 
and  its  performance  in  the  presence  of  channel  errors  is  no  worse 
than  the  nonadaptive  schemes. 

ill)  The  adaptive  classlf 1 cat  ion  scheme  is  a  variable  bit  rate  technique 
(average  bit  rate  from  one  column  to  the  next  can  vary)  which  la 
robust  (with  respect  to  channel  errors)  and  can  yield  varying  degree 
of  performances  depending  upon  the  coder  complexity.  It  ia  relatively 
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more  effective  at  low  bit  rates  (<_  1  bit/plxel)  and  its  performance 
is  degraded  as  the  image  column  size  increases. 

iv)  The  Kalman  filter  DPCM  loop  scheme  is  well  suited  for  designing 
Hybrid  coding  of  noisy  images.  As  expected,  the  performance  of  SC2 
and  separable  covariance  models  are  very  close.  The  SCI  model  per¬ 
forms  better  than  the  other  two  by  about  1  to  0.5  dB.  Three  different 
rates  have  been  considered,  viz.,  the  infinite  (no  quantizer  in  the 
DPCM  loop),  2  bits/pixel  and  1  bit/pixel  rates.  The  two  bit  rate 
quantizer  gives  a  restored  image  with  SNR  only  0.5  dB  less  than 

the  infinite  rate. 

v)  The  overall  Kalman  filter  DPCM  coder  performance  depends  on  additive 
noise  vs.  quantizer  noise  (bit  rate)  considerations.  At  low  signal 
to  noise  ratios,  the  additive  noise  dominates,  and  moderate  changes 
in  bit  rates  have  only  marginal  effect  on  improving  the  SNR  of  the 
encoded  signal.  (Note  that  the  SNR  is  bounded  by  the  Kalman  filter 
performance  without  any  quantizer  noise).  At  high  signal  to  noise 
ratios,  the  effect  of  quantization  is  more  visible. 
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'If  -  One  Dimensionol  Uni  lory  Transform 


Figure  J.l  Hybrid  Coding  Syeteo 
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Figure  3.3  Realization  of  Seralcausal  Model  Decomposition 


Figure  3.4  Adaptive  Variance  Eatlnatlon  Hybrid  Coding  Syate» 
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Figure  3.3  Adaptive  Cleeelflcatlon  Hybrid  Coding  Sjrete* 
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Figure  3.10  Classification  Map  of  the  Left  Quarter  Portion  16  *  32 
Transformed  Sub-Columns  of  an  Image  (a)  Cirl  Image 
C  (b)  Chemical  Plant  Image,  where  1  Indicates  the  Highest 

Activity  Sub-Columns,  4  Indicates  the  Lowest  Activity 
Sub-Columns . 
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Figure  3.12  CMOailrw  of  Adapstvr  Claeatf  l cat  I  on  Hybrid  Cod  I  nr.  Fenulte  for  Varloua  tauipe 
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Fleur*  3.20  Hybrid  Coding  Results  for  SCI  Model  of  the  Chemical  Plant  Image  Encoded  via 
16  «  236  Image  Scrip,  (a)  Race  Distortion  Curve  (Bit  Race  versus  Normal¬ 
ised  Kean  Square  Error)  (b)  Bit  Rate  versus  Signal  to  Noise  Ratio 


(o)  Cb) 


Figure  1.21  Hybrid  Coding  Results  for  SC 2  Model  of  the  Chemical  Plant  Image  Encoded  via 
16  *  256  lauf.*  Strip,  (a)  Rat*  Distortion  Curve  (Bit  Rat*  versus  Normalised 
Mean  Square  Error)  (b)  Bit  Rat*  versus  Signal  to  Hols*  Ratio 
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(a)  SCI  Model,  Nonadaptlve 


(d)  SC2  Model,  Nonadaptlve 


(e)  SC2  Model,  Adaptive 
Variance  Estimation 


(b)  SCI  Model,  Adaptive 
Variance  Estimation 


(c)  SCI  Model,  Adaptive  (e)  SC2  Model,  Adaptive 

Classification  Classification 

Figure  3.22  Hybrid  Encoded  Girl  Image  at  1  Bit/Pixcl  Approximately 


(a)  SCI  Model,  Nonadaptive  (d)  SC2  Model,  Nonadaptlve 

(2  Bits /Pixel)  (2  Bits/Pixel) 


(b)  SCI  Model,  Adapt ive  Variance 
Estimation  (2  Bits/Pixel) 


(e)  SC2  Model, Adaptive  Variance 
Estimation  (2  Bits/Pixel) 


(c)  SCI  Model,  Adaptive  Class!-  (f)  SC2  Model,  Adaptive  Classi¬ 
fication  (1.  7  Rlts/Pixel)  fication  (l.f>7  Rits/Pixel) 

Figure  3.23  Hvhrid  Encoded  Ctrl  Image  at  1 . 7  to  2  Rits/Pixel  Approx imatelv 
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(c)  SCI  Model,  Adaptive  (f)  SC2  Model,  Adaptive 

Cl  ass l f teat  ion  Classification 

Figure  1.24  Hybrid  Kncoded  Chemical  Plant  ImaRe  at  1  Bit/Pixel 
Approximately 
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(a)  Original  Cirl 


(d)  Original  Chemical  Plant 


(b)  Actual  Measurement  Model 
0.97  Bit/Pixel 


(e)  Actual  Measurement  Model 
1.02  Bits/Pixel 


(c)  Actual  Measurement  Model,  (f)  Actual  Measurement  Model 

2.02  Bits/Pixcl  2.06  Bits/Ptxel 

Figure  3.26  Actual  Measurement  Model  Hybrid  Encoded  Image  at  1  and  2 
Bits/Pixel  Approximately 
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Figure  3.27  Model  of  a  Binary  Symmetric  Channel 
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Figure  J.78  r.ffect  of  Channel  Error 
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(a)  SCI  Model,  Nonadaptlve 


(d)  SC2  Model,  Nonadaptlve 


(b)  SCI  Model,  Adaptive  (e)  SC2  Model,  Adaptive 

Variance  Estimation  Variance  Estimation 


Figure  3.29  Hvhrld  Encoded  Girl  Image  with  Channel  Error  P  ■  0.0001  at 
1  Blt/Plxel  Approximately 


(c)  SCI  Model,  Adaptive 
Class! f lcat Ion 


(f)  SC2  Model,  Adaptive 
Classification 


(a)  SCI  Model,  Nonadaptlve 


(d)  SC2  Model,  Nonadaptive 


(e)  SC2  Model,  Adaptive 
Variance  Estimation 


(b)  SCI  Model,  Adaptive 
Variance  Estimation 
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(c)  SCI  Model,  Adaptive 
Class! f leat Ion 


Cf)  SC2  Model,  Adaptive 
Classl f leat ion 


0.001  at 


Figure  3.30  Hybrid  F.nroded  Clrl  Image  with  Channel  F.rror  P 
1  Bit /Pixel  Approximately  1 


Flgur*  J.J2  Kilun  Fllt«r  DPC.M  Loop  for  th«  1-th  Channel  In  Hybrid  Coding 


Figure  3.33  Restored  Images  of  the  Kalman  Filter  DPCM  Loop  Scheme 
Implemented  on  a  Noisy  Girl  Image  (S/N  “  5) 


(c)  SCI  Model,  1  Bit/Pixel 


(f)  SC 2  Model,  1  Blt/Plxel 
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TaSla  1.2  III  All-Nation  F4ti«n.i  on  Equti  lu.i  (1.2.2-21),  Ctrl  liaq* 


0>)  SCI  Modal 


ofcara  ■  I  l»*l  kit  *1  Watte* 

(aft  lntlMin  tkn  w»r»»l  Ittepi  a»  n 

a*t  kli  calnltlH  ala  hW  kit  alleratlm  altorlth* 


Table  3,3  Parameters  of  Image  Models  Used  in  Hybrid  Coding 


Separable  Covariance  Model 

P.  ■  p  “  0.95 
n  v 

S2(t>  -  (l-pJl/Xj 

xi  -  l¥BV'l11 
Semlcausal  Model  SCI 

a  -  .4275,  y  -  0.1415,  B2  -  0.0198  (Girl  Image) 

a  *■  .4275,  y  «  0.1420,  p2  ■  0.0169  (Chemical  Plant  Image) 

X  ■  1  -  2a  cos  ,  1  <  i  <  N 

ci  N  —  — 

2  2  2 
B  (i)  -  BVX  « 

Pi  •  TAcl 

Semlcausal  Model _SC2 

Ph  -  Pv  -  0.95,  -‘-Pv/(1+P^> 

X  -  1  -  2a  cos  ^r-}~  -  ,  1  <  i  <  N 

ci  N  —  — 

B2  -  (1-p2)(1-p2)/(1+p2) 

n  V  /  V 

82<0  -  S2/icl 

P1  '  ph 

Actual  Measurement  Model 


Vj  (i)  -  P^  jd)  +  e.j(i) 


pi  ’  *yi)Vi<1>  *  BS  j2  XVl<1)VJ-i(1> 

2  9  99  (1-P?)  M  16  .  „ 

eUi)  -  F.e2(D  -  (l-p‘)Ev2(i)  =  -Jfij-  J^tvjd)!2 

k 

where  v^(i)  “  ith  element  of  the  jth  column  of  the  kth  16  x  M  site  image 
atrip  (k*l,...,16  here). 
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Table  1.4  Hyhtld  Coding  Results  of  (hr  Cltl  Image  Encoded  via  16  "  246  Image  Strip 
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Table  3.9  Discrete  Kalman  Kilter  and  DPCM  Algorithms 
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NONCAUSAL  MODELS  AND  TRANSFORM  CODING 

4.1  Introduction 

Transform  coding  has  attracted  considerable  attention  In  the 
recent  years  due  to  its  applications  in  image  storage  and  transmission. 
This  scheme  is  different  from  the  classical  forms  of  image  coding  where 
either  the  image  samples  are  directly  coded  (PCM)  or  the  result  of  a 
local  operation  on  pixels  is  coded  (DPCM  and  predictive  coding)  [  9) . 

In  transform  coding  the  image  is  first  unitary  transformed  and  then  the 
transformed  samples  arc  quantized  and  coded  for  transmission  or  storage. 
Fig.  4.1  shows  a  typical  image  transform  coding  system. 

Transform  coding  has  some  special  characteristics,  viz., 

(a)  Energy  Packing:  The  energy  of  samples  in  the  transform  donain 
is  generally  concentrated  near  the  low  spatial  frequences.  This  property 
is  exploited  to  achieve  data  compression. 

(b)  Channel  Error  Immunity;  Transform  coded  images  are  robust  with 
respect  to  channel  errors.  This  results  from  the  inherent  averaging  oper¬ 
ation  of  the  transform.  F.ach  reconstructed  image  sample  is  a  weighted 
function  of  all  transform  samples.  Hence,  the  channel  error  is  distri¬ 
buted  over  all  of  the  reconstructed  image  samples  so  that  it  is  less 
objectionable  to  a  human  viewer. 

In  this  chapter,  noncausal  image  representations  discussed  in 
Chapter  2  are  shown  to  yield  transform  coding  algorithms.  Adaptive 
transform  coding  scheme  analogous  to  the  classification  Hybrid  coding 
method  can  be  designed  to  account  for  changes  in  image  statistics. 
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A. 2  Image  Transform 


107 


and  the  Parseval's  energy  conservation  rule 


N  N 


N  N 


l  l  h  J  -  1  l  Ivk  £| 

i-1  J-l  1,3  k-1  £-1  K’1' 


(4.2-7) 


is  satisfied.  Common  examples  of  transform  matrices  are  the  discrete 

Fourier  (F),  Sine  (V  ),  Cosine  (Y  ),  Hadamard  (H) ,  Slant  (SL) ,  Karhunen 

s  c 

Loeve  (Y),  etc.  These  are  defined  as  follows: 


a.  Fourier  Transform  (  3  ).  The  Fourier  transform  matrix  is  defined 


as 


1  J-2iri(k-l)U-l)  \ 

t.i  "  /R  exp1 - n - / 


1  <  k,l  <  N  (4.2-8) 


where  i  *  S-l ,  and  the  discrete  Fourier  transform  of  an  image  is 
given  by 

N  N  f  \ 

v  “  k  1  1  rexP  s --==•(  (m-1)  (k-1)  +  (n-1) (1-1) ]f .  (4.2-9a) 

™,n  N  k-1  1-1  N  7 


Its  inverse  transform  is 


\,t  *  N  l  ♦  (n-l)(Ul)]| 

m-1  n-1  '  ' 


(4 . 2-9b) 


b.  Sine  Transform  (  27  1.  The  Sine  transform  matrix  Yr  is  symmetric  and 


is  defined  as 


-  r^ 

1  v  w+1 


Sin^y  1  <  k,t  <  N  . 


(4.2-10) 


It  has  been  shown  (  30  ],  if  Q  is  an  arbitrary,  symmetric,  tridia¬ 
gonal,  Toeplitt  matrix 
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(4.2-10a) 


then  y  contains  the  eigenvectors  of  Q  and  satisfies  the  relation 
s  s 


q  y  -  y  a 

8  S  8  S 


(4.2-10b) 


where  A  is  a  diagonal  matrix  whose  entries  are  the  eigenvalues  of 


Qs,  given  by 


X  *  a  +  2b  cos~-  1  <  i  <  N 
N+l  —  — 


(4. 2- 10c) 


Cosine  Transform  [  1  ].  The  Cosine  transform  matrix  y  is  defined  as 

c 


k  -  1,  1  <  fc  <  N 


(4.2-11) 


yj  ,  2£k<.v.  !<£<». 


It  has  been  shown  [  29  ] ,  the  Cosine  transform  matrix  contains  the 
eigenvectors  of  an  arbitrary  symmetric  tridiagonal  matrix  Qc  of 
order  N. 

r  a-b  b 


b  a 


(4.2-lla) 


a  b 


b  a-b 


and  satisfies  the  relation 


q  y  -  y  a 

'c  r  r  r 


(4.2-llb) 


where  Ap  is  a  diagonal  matrix  whose  entries  are  given  by 


Xci  ”  a  +  2b  cos  1  i  1  5  N  •  (4.2-llc) 

Hadamard  Transform  [54  ).  The  Hadamard  transform  matrix  is  a 
square  array  whose  entry  takes  only  binary  values  *1  and  rows  and 
columns  are  orthogonal.  The  lowest  size  2*2  orthogonal  Hadamard 
matrix  is  given  by 


(4.2-12) 


The  construction  rule  for  a  Hadamard  matrix  of  size  N  =  2n, 
where  n  Is  an  integer,  is  simple.  Let  H^,  be  a  Hadamard  matrix  of 
size  N,  the  matrix 

H2N-7r 

is  a  Hadamard  matrix  of  size  2N. 


Slant  Transform  [  55  ].  For  matrix  size  of  N  -  2,  the  Slant  trans¬ 
form  matrix  is  identical  to  the  Hadamard  transform  of  order  two. 

That  is, 


(4.2-14) 


The  Slant  matrix  of  orderN  (N»2n,n"3,4,5,* • * )  can  be  generated  in 
terms  of  the  Slant  matrix  of  order  by  the  recursive  relation 


» 


where  11^  represents  a  k  *  k  identity  matrix.  The  constants  and 
bj,  may  be  computed  from  the  recursive  relation 


bN  -  [1  +  4(aN/2>2l~li  (4.2-16) 

aN  "  2bNaN/2  ’ 


f.  Karhunen  Loeve  Transform  (  37,42  ].  Here,  we  consider  a  special 

transform  which  depends  on  the  image  statistics. 

If  (x^,  i**l,...,N)  is  a  one-dimensional  first-order  Markov 
process  with  covariance  matrix  R,  where 
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Then  the  elements  of  the  KL  transform  matrix  ¥  are  given  by  [  56  ) 
(N  ■  even) 


si 


(4.2-18.,) 


where 


2  1  -  p* 

v  •  — "  — —  ■*  ■  ■ — — — 

i  1  -  2p  cos  + 


(4. 2- 18b) 


and  (o)j)  are  the  positive  roots  of  the  transcendental  equation 


tan  n  to 


(1  -  p7)  si n _n_w _ 

cos  co  -  2p  +  p7^  cos  co 


(4.2-18c) 


Due  to  the  non-periodicity  of  the  sine  terms  in  (4.2-18a),  there  is 
no  fast  algorithm  available  for  implementing  the  Kl.  transform  matrix. 
In  matrix  form,  ¥  satisfies 

R¥  -  ¥  T  (4.2-19) 

where  T  ■  {y^}  is  the  diagonal  matrix  of  eigenvalues  of  R,  which  can 
be  calculated  via  (4.2-18c). 


4.3  Stochastic  Decoupling  of  Noncausal  Models  by  Sine  Transform  [  31  ) 

In  Chapter  2  we  discussed  three  different  noncausal  models.  Their 
spatial  structures  are  depicted  in  Fig.  4.2  together  with  the  boundary 
elements  needed  in  representation  of  an  N  x  N  square  image.  Now  we  show 
that  each  model  yields  a  stochastic  decomposition  by  the  Sine  transform. 


4.3.1  NCI  Model 


From  Chapter  2,  we  can  write  the  NCI  model  (2.5. 3-1)  in  matrix  form. 


as 


UQ+QU-c  +  B 

8  S 


(4. 3.1-1) 


where  U  and  e  are  N  *  N  matrices,  and  is  a  tridiagonal,  symmetric, 
Toeplitz  matrix  defined  in  (4.2-10a)  with  a  -  >$,  b  «  -a.  From  (4.2-10c), 
its  eigenvalues  are 


XSl  "  H  "  2“  COS  j£l. 


1  <  i  <  N  . 


(4. 3. 1-2) 


The  matrix  B  is  a  sparse  matrix  of  the  form 


B  -  Ct 


0 ,  N- 1 


U0,N+U1,N+1  ^ 
U2 , N+l 
U3,N+1 


Vi.o  ; _ _ _  J  Vl,N+l 

[_  “N+l  ,1+UN,0  “n+1,2  UN+1 , 3  •  •  •  Vl,N-l  Vl  .N^.N+l  J 


(4. 3. 1-3) 


Note  that  the  four  comer  boundary  elements  u^,  uQ  N+^,  Q  and 

i»n+i  jj+i  «re  not  needed  in  B  (see  Fig.  4.2).  Alternatively  we  write 
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(4 . 3.1-4) 


B  e  Sj  (u)  . 

—  2 

Let  us  denote  A  as  an  N  *  1  vector  of  lexicographic  ordered  elements 
of  an  N  *  N  matrix  A.  Then  (4. 3. 1-1)  has  a  Kronecker  product  form 

JU  -  (I  ©  Q  +Q„  ©OU-C  +  B  .  (4. 3.1-5) 

8  S 

From  (2.5. 3-2),  the  covariance  matrix  of  {e  }  (see  Table  2.1  under 

1  *  J 

NCI)  can  be  written  as 

R  -  EEC1  -B2(I©Q  +  Q  ©  I)  -  B2J.  (4. 3.1-6) 

s  ^  s  ^  ^ 

where  Q  is  defined  in  (4.2-10a)  with  a  -  *s,  b  “  -a  .  Now  from  (4. 3. 1-5) 
81  1 
we  obtain  the  decomposition 

U  -  U°  +  U  or  U  -  U°  +  Ub  (4. 3. 1-7) 

where  U°  -  J_1C  ,  Uv  =  J_1B  .  (4.3.1-7a) 

b 

It  is  clear  that  U,  is  obtained  by  passing  the  boundary  elements  B 
b 

through  a  noncausal  filter  J  Therefore  (4.3.l~7)  shows  that  a  dis¬ 
crete  random  field  defined  by  a  noncausal  representation  has  a  stochastic 
decomposition  in  which  Ub  is  a  random  field  completely  determined  by 

Let  denote  the  space  of  all  linear  boundary  elements  of  a  random 
field  (x.  .)  which  is  described  as 

*  •  J 

Sk(x)  -  ^l*1.k,J,Vkl]'xi,l-k,Xi,N+k1  :  1_k  1  N+k} 

where  JT  denotes  any  linear  operator. 

114 


the  boundary  elements.  Fig.  4.3  shows  the  realization  algorithm  for 
noncausal  node!  decomposition  of  (4. 3. 1-7). 


The  covariance  matrix  of  U  is  obtained  as 
X 

Ro  ^  E  U°U°  -  j'iE[c  tT]j'1  -  j'1KcJ_1  -  B2J-1J1J_1  (4. 3. 1-8) 

-  02 1 1  *  Q  +  Q  ®  I  ]-1  [  I  ©  Q  +  Q  ®  I][l  ©Qb  +  Qr©  I]"1. 

. r  fl  ^  o  ^  h  o 

If  we  denote  Y  to  be  the  Sine  transform  matrix,  then  from  (4.2-10)  it 
s 

diagonalizes  Q  .  It  follows  that  V  should  also  reduce  J,  J_,  R  ,  etc. 

s  s  1  o 

all  to  their  diagonal  form.  In  fact,  it  is  trivial  to  check  that 


n  -  (V  0  V  )  X (V  ©  Y  )  ,  for  X  -  J,J1tR.R_  (4. 3. 1-9) 

S  b  oh  1  O  L 

2  2 

where  D  is  an  N  *  N  diagonal  matrix. 

Hence  (Y  0Y  )  is  the  KL  transform  of  U°,  i.e., 
s  s 

V°  -  (V  ©Y  )U°  or  V°  -  Y  U°Y  .  (4.3.1-10) 

s  s  s  s 

Then  V°  is  the  Kl.  transformed  matrix  and  its  elements  are  uncorrelated, 
i.e., 


E  V°  V° 

i.lk.l 


"  di.)6i.V6J,t 


1  <  i,j  <  N 


(4.3.1-11) 


Therefore,  it  is  obvious  that  one  can  write  [  7  ) 

D  •  B2(  I  ©  A  +  A  ®  I)-1  (I  ®  A  +  A  ©  I)  (I  ©  A  +  A  ©I)'1 

8  8  R  j  S  j  8  8 

(4.3.1-12) 

which  gives 
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1  £  i.j  <  N 


(4.3.1-13) 


where  Xg  is  obtained  from  (4. 3.1-2)  with  a  replaced  by  . 

For  a  given  image  the  elements  V°  .  can  be  easily  realized  via  the 

MJ 


equation 


i.j  X  +  X 

Si  Sj 


(4.3.1-14) 


where 


e  -  y  c  y  ,  v°  -  y  u°y  . 

s  s  s  s 


(4.3. l-14a) 


Using  the  representation  (2.5. 3-1)  and  boundary  matrix  (B  .},  find  the 

i  t  J 

matrix  (c  .),  then  from  (4. 3.1-2)  find  the  values  of  X  and  X  ,  and 
i.j  Bj  Sj 

compute  the  values  of  e.  via  e  -  V  C  f  .  Since  y  is  the  Sine  trans- 

i.j  s  s  s 

form,  a  fast  algorithm  can  be  used  in  computing  all  the  equations  dis¬ 
cussed  in  the  above. 

4.3.2  NC2  Model 

be  written  in  matrix  form  as 

-  C  +  B  (4. 3.2-1) 

i  -  J 

I1-Jl  "  1  (4. 3.2-2) 

|i-j|  -  2 

otherwise 

and  Qg  is  a  tridiagonal,  symmetric,  Toeplitz  matrix  defined  in  (4.2-10a) 
with  a  ■  %,  b  *  -a,  and  B  is  an  N  *  N  sparse  matrix  which  contains  the 
boundary  elements  and  is  defined  in  the  following  equation. 


NC2  representat ion  (2 . 5 . 3-6)  can 


UH  +  HU+2QUQ 
s  s 


where 


i.j 


\  +  2cT  , 


-a 

a2 

0 
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An  easier  notation  for  B  is  B  c  SjCu)  (see  page  114  footnote). 

In  Kronecker  product  notation,  (4. 3. 2-1)  becomes 

(I®H  +  H©I+2Qs0QslU-e  +  B  (4. 3.2-4) 

J  U -  C  +B  (4 . 3. 2-4a) 

where  J  «  (I©H  +  H®  I  +  2Q  ®Q  ]  and  the  covariance  matrix  of  c  is 

S  8 

Rf  -  F.  ccT  -  (4. 3.2-5) 

where  a  is  replaced  by  ci^  in  and  F.qn.  (4. 3.2-4)  can  be  rewritten  as 

U  -  U°  +  Ub  (4. 3. 2-6) 

where  U°  -  j'1?  ,  UL  -  J_1B  (4.3.2-6a) 

l> 

and  the  covariance  matrix  of  tf°  is 

R0  -  F.  U°U°T  -  J_1E(€CT]  J_1  -  B2J"1J1J_1  .  (4. 3. 2-7) 

Now,  the  matrix  R0  cannot  be  diagonalized  by  the  Sine  transform  matrix 

iji  (see  (4.2-10))  and  the  Kl.  transform  of  l’°  has  to  be  the  biorthogonal 
8 

transform  which  simultaneously  reduces  H  and  matrices  to  their  dia¬ 
gonal  form. 

It  is  noted  that  H  can  be  further  decomposed  as 

H  -  Q*  ♦  Hb  (4. 3.2-8) 

where  Hb  is  a  sparse  matrix  of  only  two  nonzero  elements 


8 


(4. 3.2-9) 


Recalling  the  requirement  of  jaj  <  ^  for  NC2  model  (see  (2.5. 3-6)), 
we  nay  assume  s  [0]  and  a  reasonable  approximation  of  H  is 


H  -  Q  . 
s 


(4.3.2-10) 


Substituting  (4.3.2-10)  into  (4.3.2-14a)  yields 


J-I0H  +  H0I  +  2Q  ®  Q 

fi  s 


■  1©(T  +  Q~©I  +  2Q  ©  Q 

8  n  S  o 


(4.3.2-11) 


If  Y  denotes  the  Sine  transform,  it  is  trivial  to  prove  that 
s 

(Y  ©  Y  >  J'  (Y  ©  V  )  -  1  ©  A2  +  A2  ©  I  +  2A  ©  A  (4.3.2-12) 

ss  ss  ss  ss 

where  A  is  an  N  *  N  diagonal  matrix  with  entries  defined  in  (4. 3. 2-2). 
s 

Following  the  Bane  steps  as  we  developed  in  the  case  of  NCI  model, 
the  Sine  transform  is  a  good  approximation  to  the  biorthogonal  trans¬ 


formation  in  the  case  of  NC.2  model. 


E  V?.  S  B2 


[v 

fx  +i  l4 

L  8jJ 


(4.3.2-13) 
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Furthermore,  since  is  independent  of  a,  a  fast  Sine  transform 


exists  and  the  values  of  V  can  be  easily  realized. 

*»  J 


4.3.3  NC3  Model 

NC3  representation  is  shown  in  (2.5.3-10).  In  matrix  notation,  it 
becomes 

Q.UQ  ■  E  +  B  (4. 3. 3-1) 

»  s 


where  Qg  is  a  tridiagonal,  symmetric,  Toeplitz  matrix  defined  in 
(4.2-10a)  with  a  -  1,  b  ■  -a,  and  eigenvalues  Ag  as  defined  in  (4. 3.1-2) 

and  B  -  +  ciB2Qg  -  a2B3  (4. 3. 3-2) 


where 


U01 

U02 

c 

o 

z 

B1  - 

o 

_Vi,i 

Vl.2 

•••  Vl.N 

uio ; 

i 

!  U1.N+1 

B2  ’ 

u20  ' 

1 

l 

e 

1 

•  i 

o 

1  U2,N+1 

•  t 

“n.o  ' 

I 

1  '‘n.n+i 

r>j 

1_Vn+_i 

U0,N+1 


l - 

J  Vn.N+l 


(4.3.3-  3a) 


(4.3.3-3b) 


(4. 3. 3- 3c) 
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Thus 


B  e  Sj(u) 


(4. 3. 3-4) 


where  the  definition  of  see  page  114  footnote. 

Rewriting  (4. 3. 3-1)  in  Kronecker  product  form,  we  have 

(lj  ®Q  )U  -  JU-  £  +  B  (4. 3. 3-5) 

8  8 

where  J  ■  Q  0  Q  (4. 3. 2-6) 

8  S 

so  that 

U  -  U°  +  V.  (4. 3. 3-7) 

D 

where  U°  -  J~lc  ,  Uv  ^  J_1B  .  (4.3.3-7a) 

b 

From  Table  2-1,  it  is  noted  that  the  covariance  matrix  of  C  can  be  written 
as 

Rj.  -  E  ccT  -  P2(Qb©Qb)  •  (4. 3. 3-8) 

Thus  the  covariance  matrix  of  U11  is 

RQ  -  E  U°U°T  -  j"1E[GCT)j"1  -  B2J_1J  J"1  -  B2J_1  .  (4.3. 3-9) 

Let  be  the  Sine  transform  independent  of  a,  recalling  the  mathematical 
development  described  for  NCI  and  NC2  models,  it  is  seen  that  (Yg  ®  Y^) 
is  the  KL  transform  of  U°.  Let  us  define 


(4.3.3-10) 


(4.3.3-11) 
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Thus  the  above  equation  is  in  a  form  which  is  easily  realizable. 


For  a  given  image,  the  matrix  (B  )  is  first  computed  and  leads  to  the 

*»  J 

matrix  (U  }.  Then  from  (4. 3.1-2)  the  values  of  X.  and  A,  are  found. 

J  l  j 

A  two-dimensional  fast  Sine  transform  is  performed  over  (u°  ).  (4.3.3-11) 

i»  J 

is  used  to  evaluate  the  variances  of  each  point  of  the  transformed  image 


fv°  ) 

1  i.r  * 


4.4  Stochastic  Formulation  of  Noncausal  Models  by  Cosine  Transform 
The  Sine  transform  could  be  used  to  analyze  the  noncausal  repre¬ 
sentations  as  described  in  the  previous  section.  For  certain  special 
conditions  on  the  outmost  boundary  variables  (refer  to  Fig.  4.2),  the 
resulting  representations  could  also  be  analyzed  by  the  Cosine  transform. 

For  many  image  fields,  the  adjacent  elements  of  a  column  (or  a  row) 
are  highly  correlated  and  often  the  correlation  parameter  has  a  value 
close  to  unity.  Let  us  first  consider  the  NCI  model. 

For  a  finite  N  x  N  image  random  field,  the  boundary  elements  of  NCI 
model  that  need  to  be  specified  in  (2. 3. 3-1)  satisfy  the  conditions 


u  .  ■  u,  , 
o.J  i,J 


Vl.j  ”  UN, J 


ui.o"  “i.l 


Ui,N+l  "  Ui,N 


1  <  J  £  N 
1  <  j  <  N 
1  <  i  <  N 
1  <  J  <  N 


(4.4-1) 


This  means  the  two  outmost  boundary  rows  of  four  edges  of  the 
(N+2)  *  (N+2)  field  are  equal.  The  NCI  image  field  subject  to  (4.4-1)  is 
no  longer  a  stationary  field.  Rewriting  (2. 5. 3-1)  in  matrix  form  with 
the  constraint  (4.4-1),  we  have 


Q  U  +  U  Q  -  c 

c  c 
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(4.4-2) 


and  eigenvalues  of  Q c  are  Riven  by 

X  "  *4-20  cos  1  <  i  <  N  .  (4.4-3a) 

C  N 

The  subscript  c  indicates  Qc  is  now  diagonalized  by  the  Cosine 
transform.  In  Kronecker  product  form  (4.4-2)  becomes 


(I  ©  Qc  +  Qc  ©  I)U  "  C  .  (4.4-4) 

Therefore  U  •  J  (4.4-5) 

where  J"(I©QC  +  QC©I)«  (4. 4- 5a) 

The  covariance  matrix  of  R  is 

R  -  E  UUT  -  J-1E[ecT]  j"1  -  J-1R  j'1  (4.4-6) 

where  R  •  E  E ?  .  (4.4-6a) 


Eqn.  (4.4-5)  does  not  involve  a  term  in  U^,  compare  with  (4. 3. 1-7),  any 
more. 

Since  we  have  used  the  boundary  conditions  of  (4.4-1)  in  deriving 
(4.4-2),  it  means  that  the  error  signal  (e^)  along  the  four  boundary 
edges  is  no  longer  a  completely  arbitrary  random  variable.  For  stationary 
NCI  model,  the  covariance  matrix  of  C  was  defined  in  (4. 3. 1-6),  where 
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(4.4-8) 


The  above  equation  could  be  diagonalized  by  Cosine  transform,  and 
the  following  expression  Is  obtained  for  Rf 

R_  -  82(Q  ©  I  +  1  ©  Q  )  -  P2J.  (4.4-9) 

1  C1  1 

where  J.  -  (Q  ©I  +  I®Q  )  .  (4.4.9a) 

l  Cj 

Substituting  (4.4-9)  to  (4.4-6),  yields 

R  -  .  (4.4-10) 

Following  the  same  procedure  as  we  did  In  section  4.3.1,  It  is  easy  to 

show  (Y  ©  ¥  ), where  Y  Is  the  Cosine  transform  matrix,  could  reduce  J 
c  c  c 

and  Jj  to  their  diagonal  forms.  Therefore 
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mm. 


E  VVT  '■  e2(A„0  I  +  I©A,)_1(A  0  I  +  I  ®  A  XA^  ®  I  +  I  0  A  )~1 


C1  C 


which  implies 


2  „  2  (  \  XClj 

Vi;j " 8  Tx  7~Ty  1£1’j£N 

l  Ci  Cj / 


(4.4-11) 


(4.4-12) 


where  (X  }  and  (X  }  could  be  computed  from  (4. 4- 3a).  It  is  obvious 
CH  ci 

that  the  elements  v^  ^  can  be  quite  easily  realized  and  given  by 


'i.j  x  +  x 

ci  cj 


(4.4-13) 


where 


e  -  ?  e  ?  ,  V  S  T  ITT 

c  c  c  c 


(4.4-13a) 


Once  again,  fast  algorithm  can  be  used  to  implement  the  Cosine 
transform  in  the  above  equations. 

For  NC2  and  NC3,  similar  models  can  be  obtained  by  simply  assuming 
the  outermost  boundary  values  to  be  equal.  Table  4.1  summarizes  the 
formulas  for  the  three  models  which  are  useful  in  the  design  of  transform 
coders  for  random  fields  represented  by  noncausal  models. 


4.5  Transform  Coder  Design  for  Noncausal  Models  via  Cosine  Transform 
In  coding  problems,  one  would  like  to  find  the  minimum  number  of 
bits  required  to  represent  a  random  field  subject  to  the  specific  amount 
of  distortion.  A  common  approach  of  evaluating  the  performance  of  a 
coding  scheme  is  to  obtain  its  distortion  vs.  rate  characteristics. 
Referring  to  Fig.  4.1,  which  shows  the  overall  transform  coding  system, 
an  image  U  is  unitarlly  transformed  to  obtain  V.  Each  element  of  V  is 
then  quantized  Independently  by  a  zero  memory  quantizer  and  is  used  for 
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transmission  and/or  storage.  For  reconstruction  of  the  original  image, 
one  simply  performs  the  inverse  transformation  of  the  quantized  samples. 

One  basic  property  of  the  KL  transform  is  that  it  completely  de¬ 
correlates  the  transform  samples.  Other  unitary  transforms  reduce  the 
intei^sample  correlation.  The  transform  coder  efficiency  is  maximized 
when  the  transformed  samples  are  uncorrelated.  The  quantizer  design 
depends  on  the  probability  distribution  of  the  transformed  samples.  In 


our  study,  we  have  found  the  transform  sample 
well  by  the  Gaussian  density  model 


v 


i  *  j 


can  be  modeled  quite 


p(x)  -  exp[-(x-y)2/2o2) 


(4.5-1) 


where  u  and  0  are  the  mean  and  standard  deviation  of  the  random  variable 
x.  Many  times  the  first  sample  ^  is  modeled  by  a  Rayleigh  density, 
because,  for  many  transforms,  it  is  proportional  to  the  average  value  (or 
the  so  called  d.c)  of  the  image  data,  which  for  non-negative  image  data 
would  be  a  non-negative  random  variable.  However,  if  the  image  data  has 
been  converted  to  a  zero  mean  data  before  processing,  then  the  Gaussian 
density  model  suffices.  In  the  sequel,  without  loss  of  generality,  we 
will  assume  (u.  ,)  to  be  zero  mean  random  variables,  therefore,  {v  ) 

i.J  ro.n 

arc  also  zero  mean  random  variables.  Let 


02  -  E  v2  (4.5-2) 

v  m,n 

m,n 

be  the  variance  of  the  transformed  sample  v  .  In  Fig.  4.1  if  there 

m,n 

are  no  channel  errors,  then  the  coder/decoder  will  generate  a  perfect 

transaction  and  we  will  have  v*  -  v*  The  average  mean  square  dis- 

rn  |  n  oi  fH 

tortion  in  encoding  entire  N  *  N  image  is  defined  as 


N  N 


D  mJ2  l  l  E  (G*i  j"ul  1>2‘ 
n  i-1  j-1  ,:1  1,3 


(4.5-3) 

Since  the  transformation  is  unitary,  the  Parseval's  relation  implies 


1  S  - 

D  -  -r  £  l  E  v  -v  ) 
„2  **,  **,  m,n  m,n 

N  m-1  n-1 


and  in  the  absence  of  channel  errors 

.  N  N 

D  -  — 2  l  l  E  (v*  -v  )‘ 
N2  m-1  n-1  m,n  “,n 


(4.5-4) 


(4.5-4) 


Assuming  all  the  quantizers  are  identical  in  their  characteristics, 
f(b),  where 

f(b)  -  mean  square  distortion  of  a  quantizer  with  2^  levels 
for  unit  variance  input  random  variable 

p  -  desired  average  bit  rate  in  bits/pixel  of  the  image  random 
field 

we  can  write 

2f(b“-n)0v(m,n)  <4*5'5) 

N  m-1  n»l 

where  b  -  number  of  bits  allocated  to  the  sample  v  .  Now,  since 
m,n  m,n 

the  total  number  of  bits  available  is  fixed,  i.e.. 


N  N 

I  l  b  -  N2p 
m-1  n-1  m,n 


(4.5-6) 


and 


b  >0 
m,n  — 


(4.5-7) 


the  overall  distortion  given  by  (4.5-5)  is  minimized  by  finding  the  opti¬ 
mal  bit  allocation  among  the  various  samples  subject  to  the  constraints 
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of  (4.5-6)  and  (4.5-7).  Following  Chapter  3,  section  3.3.2,  and  Appendix 

A,  the  bit  allocation  problem  can  be  solved  as  before  and  the  distortion 

vs.  rate  characteristics  can  be  determined.  For  the  noncausal  models  con- 
2 

sidered  here,  Ov(m,n)  are  obtained  from  Table  4.1  for  Cosine  transform 
coding.  The  overall  coding  system  is  shown  in  Fig.  4.4(a).  In  general 
when  the  conditions  of  (4.4-1)  are  invalid,  boundary  variables  are 
arbitrary,  we  can  use  the  following  method. 

4.6  Fast  KL  (Sine)  Transform  Data  Compression 

In  either  noncausal  model  (NCI,  NC2  or  NC3) ,  the  decomposition 
U  ■  U°  +  plays  an  important  role  in  the  analysis.  Fig.  4.4(b)  illus¬ 
trates  a  Sine  transform  data  compression  scheme.  First  the  boundary 

* 

samples  are  quantized  and  the  boundary  response  is  computed.  Sub¬ 
tracting  U*  from  U  gives  U°,  which  is  nearly  equal  to  U°,  if  quantization 
errors  in  boundary  samples  are  small.  The  covariance  matrix  of  11°  is  a 
function  of  Q  for  NCI  and  NC3  models  and  a  function  of  Q  and  H  for  NC2 
model.  It  has  been  shown  that  Sine  transform  serves  as  the  exact  Kl.T  for 
NCI  and  NC3  models  and  a  good  approximation  to  the  NC2  model.  U  is 
therefore  quantized  with  a  bank  of  noncausal  model  influenced  quantizers 
and  transmitted  over  the  channel. 

In  actual  implementation,  a  large  image  would  be  first  divided  into 
small  blocks  and  each  block  is  coded  independently.  Since  the  noncausal 
models  require  separate  coding  of  image  boundaries,  the  nunfrer  of  boundary 
samples  that  are  needed  to  be  quantized  for  each  block  can  be  minimized 
by  overlapping  the  boundaries  of  the  various  image  blocks,  namely,  only 
the  right  vertical  and  lower  horizontal  need  to  be  quantized.  Rate 


distortion  calculations  for  one-dimensional  noncausal  models  (with  one 
overlapping  boundary  point)  have  shown  to  give  lower  bit  rates  than  the 
conventional  Kl.T  or  other  unitary  transforn  coding  without  any  overlap. 
Details  of  these  aspects  and  comparisons  with  conventional  transform 
coding  techniques  are  given  in  [  36  ] . 

4.7  Experimental  Results 

Several  computer  experiments  have  been  performed  to  simulate  the 
transform  coding  schemes  discussed  in  the  previous  sections.  In  all  the 
experiments  only  Cosine  transform  based  scheme  was  implemented.  A 
236  x  256  image  was  divided  into  16  image  blocks,  each  of  size  64  x  64. 
Each  block  was  transform  coded  independently.  The  size  of  image  block 
could  be  chosen  arbitrarily,  the  size  64  *  64  being  chosen  for  comparison 
with  the  methods  of  the  next  chapter. 

Figure  4.5  shows  a  typical  bit  allocat ion  pattern  of  a  NCI  Cosine 
transform  image.  Tables  4.2  and  4.3  summarize  the  results  for  the  Girl 
and  Chemical  Plant  images  respectively.  Figure  4.6  shows  rate  distortion 
curve  and  rate  vs.  SNR  for  the  Girl  image.  The  corresponding  character¬ 
istics  for  the  Chemical  Plant  image  are  shown  in  Figure  4.7.  Figure 
4.8(a)-(c)  show  encoded  Girl  images  at  0.5  bit  /  pixel,  1  bit/pixel  and 
2  bits/pixel.  Figure  4.9(a)-(c)  show  encoded  Chemical  Plant  Images  at 
0.5  bit/pixel,  1  bit/pixel  and  2  bits/pixel.  As  expected,  the  perform¬ 
ance  of  this  model  (NCI)  is  consistent  with  the  conventional  transform 
coding  methods. 
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4.8  Classification  Transform  Coding 

We  have  mentioned  in  the  last  section  that  in  order  to  increase  the 
transmission  efficiency  and  ease  the  burden  of  buffer  storage,  the  test 
image  is  divided  into  snail  image  blocks.  Then  each  image  block  is 
transform  coded  independently.  In  the  classification  transform  coding, 
one  additional  step  is  added  to  allocate  flexible  bits  to  each  transform 
channel  quantizer.  It  says  each  image  block  is  classified  as  belonging 
to  one  of  K  predetermined  classes  according  to  its  own  activity.  The 
activity  in  the  image  block  is  measured  by  the  variance  of  that  image 
block.  For  a  given  class  of  images  the  probability  distribution  function 
of  all  the  image  block  variances  can  be  predetermined.  Note  that  this 
probability  distribution  function  would  depend  on  the  size  of  the  image 
block  used  in  transform  coding. 

The  transform  signal  in  each  transform  channel  is  processed  as  before 
except  that  the  signal  in  the  transform  domain  is  quantized  according  to 
the  classi f icat ion  of  that  image  block.  In  this  way,  the  transform  channel 
quantizer  is  adapted  to  the  changes  of  activity  in  the  image  blocks. 

Image  blocks  with  high  activity  are  assigned  more  quantizer  bits  than 
those  of  low  activity.  Therefore  the  quantizer  bits  are  employed  more 
efficiently.  This  scheme  is  somewhat  analogous  to  the  adaptive  classi¬ 
fication  hybrid  coding  described  in  Chapter  three. 

This  scheme  also  requires  the  classification  information  to  be  trans¬ 
mitted  to  the  receiver.  This  can  be  done  by  sending  an  extra  log^K  bits 
per  image  block  or  ^  log^K  bits/pixel  (where  K  stands  for  the  number  of 
classes  and  N  is  the  size  of  an  image  block).  F.xperimentally ,  this  value 
turns  out  to  be  very  small.  For  example,  in  the  case  of  K  **  4  and  M  -  16, 
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the  overhead  for  its  transmission  is  only 


lo*24 


16  *  16  128 


bit  /  pixel. 


Let  us  define  y  probability  of  the  kth  class  and  also  assume  for 
each  class  k,  the  inage  model  parameters  have  been  predetermined  and 
are  knovm.  Then  following  the  same  development  of  adaptive  classifica¬ 
tion  hybrid  coding,  we  can  write  the  expected  average  distortion  for 
each  image  block  as 


N  N  K 


D  *  l  l  l  °v  (k)f(bi  j  k>Pk 

N  i-1  j-1  k-1  Vi.j  k 


(4.8-1) 


where  b  ,  .  >  0  and  represents  the  ijth  transform  channel  quantizer  bit 
and  kth  class. 

Another  constraint  could  be  the  expected  average  bit  rate  must  be 
constant ,  i.e. , 


N  N  K 


l  l  l  I  bi  j  kPk  ’  P  * 

N  i-1  j-1  k-1  * 


(4.8-2) 


The  above  hit  allocation  problen  can  he  solved  along  the  same  lines 
in  the  previous  section. 

Figure  4.10  shows  a  typical  bit  allocation  pattern  for  two-class 
problem  with  eq.ial  prohahl  1  it ies  (pfc  -  0.5)  for  each  class  and  size  of 
the  image  block  was  64  *  64.  Figure  4.11(a)  shows  the  corresponding 
4x4  classification  nap  and  image  block  measure  variances.  Classifi¬ 
cation  map  and  image  block  measure  variances  for  the  Chemical  Plant 
image  with  name  kind  of  strategy  is  shown  in  Figure  4.12(b). 

The  scheme  requires  the  following  additional  steps  compared  to 
ordinary  transform  coding. 
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1.  Measurement  of  the  variance  of  each  image  block. 

2.  Classification  of  each  image  block  to  one  of  K  classes  based 
on  the  assumption  of  p^. 

3.  Storing  of  K  bit  allocation  tables,  one  for  each  class. 

Results  from  several  simulations  show  that  this  scheme  gives,  rela¬ 
tively,  a  large  improvement  in  SNR  for  lower  bit  rates  =0.5  bit  /  pixel 
and  relatively  smaller  image  block  sizes.  At  higher  bit  rates  (»2  bits/ 
pixel)  and  larger  image  block  sizes,  the  improvement  in  SNR  becomes  less 
significant  if  at  all,  and  is  evidenced  by  appearance  of  a  "blocking" 
effect  in  the  encoded  image.  This  is  because  as  the  image  block  size  is 
increased,  the  ambiguity  in  classifying  it  into  different  classes  in¬ 
creases  since  the  activity  in  a  large  image  block  is  no  longer  uniform. 
Generally,  the  following  conclusions  can  be  made. 

1.  Classification  transform  coding  scheme  can  improve  the  SNR  over  the 
standard  transform  coiling  (based  on  NCI  model)  for  simple  structured 
image  (like  Girl  image). 

2.  The  complexity  of  implementing  this  scheme  is  marginal.  Rit  allocation 
is  wholly  determined  by  the  model  with  a  proper  adjustment  of  the 
actual  image  block  variance. 

3.  As  the  number  of  classes  goes  up,  the  quality  of  the  reconstructed 
image  improves  if  the  size  of  image  block  is  much  less  than  the 
actual  size  of  the  test  image. 

Tables  4.2  and  4.3  list  the  results  of  classification  transform 
encoded  Girl  and  Chemical  Plant  images,  both  with  64  *  64  image  block  size 
and  two  classes.  The  rate  distortion  curves  of  this  scheme  are  shown  in 
Figures  4.6  and  4.7.  Finally,  the  encoded  images  at  0.5  bit  /pixel,  1 
bit/pixel  and  2  bits/pixel  are  shown  in  Figures  4.8  and  4.9. 
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Integer  Bit  Allocation  for  NCI  Model  Cosine  Transform  Encoded 
Girl  Image  with  Average  Bit  Rate  ■  1  Bit/Pixol 


(d)  0.5  Bit/Pixel, 
Classification 


(a)  0.5  Bit/Pixel,  Standard 


(e)  1  Bit/Pixel, 
Classification 


(b)  1  Blt/Pixol,  Standard 


(f)  2  Bits/Pixel, 
Classl f lcat ion 


(c)  2  Bits/Pixel,  Standard 


Figure  4.8  Cosine  Transform  F.ncoded  Girl  Image 


(a)  0.5  Bit/Pixel,  Standard 


(d)  0.5  Bit/Pixel, 
Classl f icat ion 


(b)  1  Bit/Pixel,  Standard 


(e)  1  Bit/Pixel, 
Classification 


(f)  2  Bits/Pixel, 
Classification 


(c)  2  Bits/Pixel ,  Standard 


Cosine  Transform  Encoded  Chemical  Plant  Image 
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Figure  4.10  Bit  Allocation  with  Clanni f icat ion  at  Bit  Rate  •  1  Blt/Plxel 
(k™2,  Girl  ImaRC  and  64  *  64  Image  Block  Site) 
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Figure  4.11  Calssi f icatlon  Map  and  Image  Block  Measurement  Variances 
with  Two  Classes  (Equal  Probabilities  ■  0.5) 

(a)  Girl  Image  (b)  Chemical  Plant  Image. 

where  2  Indicates  lower  activity  and  1  indicates  higher  activity, 

256  256 

size  of  the  classification  map  is  obtained  as  (— jj-)  x  and 

N  -  64  in  this  example. 
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Table  4.2  Cirl  Image,  Cosine  Transform  Coding  Results, 
NCI  Model  with  64  x  64  Image  Block 


Standard  Classification 

Bit  Rate 

N.M.S.E. (X)  SNR  N.M.S.E.(t)  SNR 


0.5 

0.8213 

30.396 

0.8208 

30.399 

1.0 

0.6399 

31.480 

0.5906 

31.828 

1.5 

0.4929 

32.614 

0.3956 

33.569 

2.0 

0.3075 

34.662 

0.2556 

35.465 

Table  4.3  Chemical  Plant  Image,  Cosine  Transform  Coding  Results, 
NCI  Model  with  64  x  64  Image  Block 


Standard  Classification 

Bit  Rate 


N.M.S.E. (X) 

SNR 

N.M.S.E. (2) 

SNR 

0.5 

1.4067 

27.214 

1.4711 

27.020 

1.0 

1.0443 

28.508 

1.0699 

28.403 

1.5 

0.8320 

29.495 

0.8087 

29.619 

2.0 

0.6106 

30.830 

0.6244 

30.742 
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CHAPTER  FIVE 


IMAGE  MODELS  FOR  FEATURE  CODING 


5.1  Introduction 

In  this  chapter  we  consider  an  image  coding  technique  which  preserves 
the  attractive  characteristics  of  standard  transform  coding,  i.e.,  energy 
packing  and  channel  error  immunity,  as  well  as  Improves  the  visual  appear¬ 
ance  of  an  encoded  image.  This  technique  is  called  Feature  transform 
coding.  The  basic  idea  behind  this  technique  lies  in  segmenting  an  image 
which  is  considered  to  be  a  composition  of  two  parts,  deterministic  and 
stochastic,  where  the  stochastic  part  is  dealt  with  by  known  statistical 
analysis  according  to  the  image  model  algorithms  developed  in  Chapter  4 
and  the  deterministic  part  (which  could  be  called  features)  is  obtained 
as  the  residual  after  the  stochastic  part  is  removed. 

Other  two  source  modeling  of  images  has  been  considered  earlier  by 
Schreiber  [  62  ] ,  Yan  and  Sakrison  [  74  ] .  The  present  approach  is  some¬ 
what  more  general  and  is  based  on  stochastic  models  of  images  by  partial 
difference  equations. 

In  section  5.2  we  demonstrate  the  concept  of  feature  extraction  via 
image  model  operators.  In  section  5.3  we  then  relate  the  transform  coding 
scheme  with  the  feature  extraction  concept  and  a  set  of  operations  per¬ 
formed  with  this  new  technique  is  described  in  detail.  Finally,  in 
section  5.4,  we  present  results  of  encoded  Images  and  compare  them  with 
standard  NG1  Cosine  transform  encoded  image  at  the  same  bit  rate. 

5.2  Feature  Extraction  via  Image  Models 

The  problem  of  feature  extraction  is  to  determine  the  rough  profile 
of  an  object  within  a  given  Image.  It  is  noted  that  most  important 


informations  of  an  image  lie  on  features  which  contribute  the  visual 
appearance  of  an  image.  Here  a  new  technique  is  presented  in  which  the 
edge  points  (feature)  is  detected  from  the  knowledge  of  image  models. 

5.2.1  Image  Modeling  Operator 

A  given  image,  U,  can  he  treated  as  a  sample  function  of  a  two- 
dimensional  random  field.  We  consider  the  random  field  as  a  composition 
of  two  parts,  the  first  part  representing  the  relatively  smooth  variations 
within  the  image  which  can  be  described  by  stochastic  image  models  and 
the  other  representing  the  unpredictable  part  or  deterministic  part  called 
the  features  or  edges.  Thus,  we  can  write  the  image  U  as 

U  -  U  +  U,  (5.2. 1-1) 

s  d 

where  L's  and  represent  the  stochastic  and  deterministic  parts  respec¬ 
tively. 

Typically,  the  stochastic  part  representing  smooth  variations  con¬ 
sists  of  background  and  low  spatial  frequency  information  and  the  deter¬ 
ministic  part  describes  the  abrupt  changes  (large  variations)  usually 
caused  by  the  distinct  edges  of  an  object.  Therefore,  a  more  effective 
coding  method  would  be  to  segment  the  image  into  the  said  parts  and  code 
them  efficiently. 

Typically,  an  infinite,  stationary  random  field  mav  be  represented 
by  its  spectral  density  function  (SDF).  Several  candidates  have  been 
proposed  in  Table  2.1.  Let  X  denote  an  operator  acting  on  a  discrete  image. 

Denoting  e  ns  the  output  signal  after  applying  this  operator,  we  can 

write 

C.  A/lU.l.  tA±Z ludl.  (5. 2. 1-2) 
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A  common  example  of  «T  Ls  the  Laplacian  operatort  which  takes  the  dif¬ 
ference  of  an  image  clement  with  an  algebraic  average  of  its  nearest 
four  neighbors.  Most  feature  extraction  (or  edge  detection)  problems 
could  be  formulated  as  a  segmentation  of  Ug  and  in  the  above  equation, 
of  course,  with  different  operators. 

Letting  t  be  the  arbitrary  selected  threshold  which  is  employed  to 
classify  e(i,j),  the  (i.j)th  component  of  G,  the  segmentation  criterion 
is 


cs(i,j) 

Cd(t.j) 


e(i.j) 

if 

|e(i,j) !  <  t 

t 

if 

| c( i , J ) 1  >  t 

0 

if 

|c(i,J)|  <  t 

c(i, J)-Sign  C(i,j)  if  |e(i,J)|  >  t 


(5. 2. 1-3) 

(5. 2. 1-4) 


Thus  we  compare  ;C(i,j)[  with  t,  and  if  |e(i,j)|  >  t,  then  the 
point  (i.j)  belongs  to  the  set  of  boundary  points;  otherwise,  it  belongs 
to  the  stationary  part.  The  values  assigned  to  Gg  and  c ^  are  according 
to  the  above  equations. 


+  The  Laplacian  of  a  cont inuous  function  f(x,y)  is  defined  as 

x  A  „2,  3zf  ,  32f 

e(x,y)  “  V  f  -  — r  +  — - 

Zx* 


and  its  discrete  equivalent  is  given  by  the  difference  equation 


Ci,J  "  (fi+l,j 


+  fi-l,J  +  fi,J+l  +  "  4fi,J 


where  f.  .  »  f(x.,y.),  e  -  e(x,y)AxAy,  Ax  and  Ay  are  the  finite 
MJ  1  J  1  *  J 

difference  intervals  in  the  x  and  y  coordinates  respectively. 
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There  are  two  simple  ways  of  selecting  the  threshold  t.  First, 
using  the  histogram  of  lc(l,j)i,  a  threshold  value  is  picked  such  that  a 
desired  percentage  of  tor.il  Image  elements  would  be  detected  as  edge 
points.  Second,  we  could  let  t  «  c 0,  c  *  constant,  O  ■  standard  devia¬ 
tion  of  |e(i,j)|.  Now,  let  us  consider  the  selection  of  X. 


(x(i,j)} 


h(i.j) 

H(2j ,z2) 

fy(i.J)) 


Let  fx(i,j)}  and  {y(i,J)}  be  the  input  and  output  sequences,  res¬ 
pectively,  of  a  two-dimensional  linear  system  with  transfer  function 
lUij.z  ).  (x(i,j)}  and  fy(i,j))  are  related  by 

Sy(zl,z2^  "  •Z2^H^Z11 ,z21^Sx^zl,z2^  *  (5. 2. 1-5) 

For  an  input  white  noise  field  (w(i,J)},  the  output  would  be  a 
stationary  random  field,  say  {s(i,))}t  such  that  its  SDF  is  given  by 

Ss(zrz2)  "  H(*1.*2)H(*j1.*21)Sw(*i.*2>  •  (5. 2. 1-6) 

For  the  sake  of  simplicity, the  SDF  of  the  white  noise  field  is  taken 
to  be  unity.  Then  (5. 3. 1-6)  becomes 

Ss(Vz2>  "  *  (5. 2. 1-7) 

Therefore  the  transfer  function  can  be  viewed  as 

H(zltz2)  -  tS8(*1,*2)]'‘  .  (5.2. 1-8) 
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Taking  a  two-dimensional  Z  transform  on  both  sides  of  (5. 2. 1-1) 
yields 

V(z1,z2)  ■  V#(z^ .Zj)  +  Vd^rl,z2^  »  (5. 2. 1-9) 

where  V  -  Z{U),  V  -  7.{l)  )  and  V.  -  Z(U.}.  (5.2.1-9a) 

S  s  d  d 

It  is  obvious  that  V  (z  ,z  )  ■  H(z, ,z  )W(z ,z,)  if  fw(i,j))  and  fu  (i,j)) 
are  the  input  and  output  of  a  two-dimensional  linear  system  with  transfer 
function  H(ZpZ  ).  Substituting  (5. 2. 1-8)  into  (5. 2. 1-9)  we  get 

V(Zj,z2)  -  [S9(z1,z2)]‘4W(z1,z2)  +Vd(zr«2)  (5.2.1-10) 

where  W(*^,z2)  is  the  two-dimensional  Z  transform  of  the  white  noise  field. 
Rearranging  the  terns  of  (5.2.1-10),  we  have 

(Ss(zj  ,z2)  pSKZj.Zj)  -  W(z j ,z2)  +  [ z j  *z2>  1 

Denot  ing  € (Zj  , ”  [Sr(z^,z2))  ,7.^)  and  taking  inverse  two-dimensional 

Z  transform  on  both  sides  of  (5.2.1-11),  we  get 

C(i,J)  -  e8(i,J)  +  Gd(i,))  V  i.j  (5.2.1-12) 

where  c  .  is  the  inverse  Z  transform  of  S  '*V.,  and  C  is  the  white  noise 
d  s  d  s 

field.  For  finite  image,  we  define  a  two-dimensional  transform  of  the 
image  as 

v£  ffUl-TUT1  (5.2.1-13) 

where  superscript  T  indicates  conjugate  and  transpose  operations  and  Y 
is  the  unitary  KI.  transform  of  U.  (The  above  equation  implies  T  is  a 
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separable  transform,  which  will  be  true  for  our  image  models).  From  the 
foregoing  discussion,  if  we  define 

JC  -  Tlsf*V  (5.2.1-14) 

where  superscript  -1  Indicates  inverse  operation  and  Sg  Is  the  SDK  of 
stationary  random  field,  then  JC  would  be  the  desired  edge  detection 
operator. 


5.2.2  Approximate  Rational  Function  Operator 
The  SDF  of  NCI  image  model  is  given  by 


Ss(zrz2)  “ 


1  -  m(z^  +  z^  + 


(1  -  afzj  +  Zj1  +  z2  +  z2*^2 


(5. 2. 2-1) 


where  ra  -  and  P  ■  1,  to  be  consistent  with  the  formula  given  in 

_L 

Table  2.1.  If  we  assume  R(z,,z.)  -  [S  (z,,z,))  then  the  alternative 

I  1  si/ 

expression  for  R(Zj,z_)  is  given  by  (see  Appendix  B) 


R(  >  I 

1  2  n-l  2"  '(n-1)! 


1  -  r  + 


r)  n-l  n 
2n|  r  X  ’ 


(5. 2. 2-2) 


By  picking  different  values  of  n,  we  obtain  various  orders  of  R(zj,z2) 
with  which  the  edge  detection  operator  could  be  defined  as 

JT  -  R(n)(zrz2)  n  -  1,2,...  (5. 2. 2-3) 

This  implies  JT  is  a  two-dimensional  nth  order  difference  operator. 

5.2.3  F.dge  Information  Subtraction  Procedure 

We  have  stated  the  concepts  of  image  modeling  operator  and  edge 
points  segmentation  criterion  in  the  preceding  sections.  The  sequence 
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of  operations  needed  to  be  performed  may  now  be  listed  as  follows: 

T 

1.  Take  discrete  unitary  transform  of  the  input  image,  V  ■  Yllf  , 
where  ^  denotes  the  chosen  discrete  unitary  transform. 

2.  Based  on  the  image  model,  find  the  proper  representation  for  the 
SDF,  S^,  in  the  V  domain. 

3.  Calculate  and  multiply  it  with  V,  yielding  €  "  SgSf. 

*1  -T 

4.  Take  inverse  discrete  unitary  transform  of  €,  obtaining  e  ”  Y  1  €  Y 

5.  Pick  a  reasonable  threshold  to  segment  f  .  from  c  according  to 

O  R 

equations  (5. 2. 1-3)  and  (5. 2. 1-4). 

6.  Generate  edge  address  map  which  is  the  set  of  (i,j)  coordinates 
of  each  detected  edge  point. 

Once  the  edge  information  is  known,  the  sta t lonnry  field  and 
deterministic  field  of  the  given  image  can  be  obtained  by  the  following 
fairly  obvious  way: 

T 

7.  Take  discrete  unitary  transform  of  C^,  getting  € ^  —  T  r  Y  . 

8.  Preraultiply  with  and  define  ■  Ss^d‘ 

j 

9.  Take  inverse  discrete  unitary'  transform  of  to  obtain  *  Y  V^Y . 

10.  Subtracting  U.  from  the  input  image  U,  gives  the  stationary  field 

a 

as  U  -  U  -  IK. 

s  a 

We  have  considered  two  transforms,  i.e.,  discrete  Fourier  transform 
and  discrete  Cosine  transform  as  the  two-dimensional  unitary  transform 
operators  in  our  preliminary  feature  extraction  study.  The  experiments 
were  performed  on  the  256  x  256  Girl  image  in  64  x  64  image  blocks.  For 
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comparison,  Che  thresholds  for  different  methods  of  edge  detection  were 
chosen  to  yield  approximately  equal  number  of  detected  edge  points.  The 


formulas  for  evaluating  the  value  (k,£)  will  be  different  according  to 
which  transform  is  actually  implemented  and  the  image  model  used.  The 
NCI  model  was  chosen  for  all  experiments.  First  we  let 


be  the  approximation  of  the  SDF  in 


Ss(k*l)  “  Sslzl  "  e 


the  discrete  Fourier  transform  domain  when  N  is  sufficiently  large.  How¬ 
ever,  in  the  case  of  Cosine  transform,  the  numerical  value  of  SB(k,t)  is 
taken  from  the  cosine  transform  variance  formula  (see  Table  4.1  ).  Fig¬ 
ures  5.1(a)  and  5.1(b)  show  the  resulting  images  of  these  two  experiments. 
A  significant  visual  difference  can  nr  observed  from  these  two  figures. 

Tlie  block  effect  obtained  with  the  discrete  Fourier  transform  approxi¬ 
mation  is  not  present  in  using  the  Cosine  transform.  Hence,  it  is 
preferable  to  use  the  Cosine  transform  in  the  study  of  feature  extraction 
concepts. 

If  the  approximate  rational  operator  k  is  used,  rather  than  , 
c(i,J)  could  be  calculated  directly  in  the  spatial  domain  (because  R^°^ 
becomes  a  two-dimensional  nth  order  difference  operator).  Segmentation 
of  e(i,J)  could  be  performed  as  before.  The  spatial  structure  of  R^n^  as 


a  difference  operator  is  shown  for  various  values  of  n  in  Appendix  B. 


5.2.4  Image  Enhancement 

A  typical  image  is  usually  made  up  of  objects  having  edges  at  which 
the  brightness  changes  abruptly.  It  has  been  known  to  investigators  in 
the  image  processing  field  that  the  human  visual  system  is  sensitive  to 
these  edges.  Bv  this  argument  a  "good"  image  should  result  if  these 
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edge  points  can  be  Isolated  and  emphasized  so  as  to  make  the  image  appear 
sharp.  We  have  initiated  a  study  towards  this  end.  Following  the  concept 
developed  in  the  previous  sections,  the  stationary  (stochastic)  field  and 
deterministic  (feature)  field  of  any  given  image  can  be  easily  generated. 
The  way  to  “sharpen"  the  image  can  be  done  by  the  following  equation 

U'  -  aUg  +  BUd  V  i,J  (5. 2.4-1) 

where  Ug  represents  the  stationary  field  which  contains  mostly  background 
Information,  and  U.  represents  the  deterministic  field  which  results  from 

U 

the  edge  information. 

Using  (5. 2. 4-1),  one  can  change  various  values  of  a  and  B  and  achieve 
different  enhancement  (sharp  or  blur)  images  U' .  Figures  5.311  and  5.411 
demonstrate  the  results.  The  thresholds  used  in  these  experiments  are 
given  in  Table  5.1.  It  is  seen  that  U g  +  l.SU^  R*ves  tbe  best  enhance¬ 
ment  image  among  the  six  combinations  in  our  study. 

The  overall  block  diagram  of  this  system  is  shown  in  Fig.  5.2,  which 
consists  of  two  functions,  decomposition  and  synthesis.  The  stochastic 
modeling  operator  is  determined  by  a  chosen  image  model,  and  (5. 2.4-1) 
is  implemented  by  adjusting  the  two  values  a  and  B. 

5.3  Transform  Coding  Using  Nonstationary  Image  Representation 

In  the  foregoing  presentation,  we  considered  a  new  concept  of  image 
decomposition,  an  image  can  be  separated  into  its  stochastic  (stationary) 
and  deterministic  (feature)  components,  i.e., 

U  «  U  +  U  .  (5.2. 1-1) 

8  a 


-! 
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One  of  the  possible  applications  of  (5. 2. 1-1)  is  to  extend  It  to 
image  coding  or  data  compression.  In  comparison  with  current  existing 
coding  schemes,  including  the  one  presented  in  Chapter  4,  this  scheme  is 
found  to  give  better  performance  in  terms  of  resolution  for  the  same  bit 
rate,  and  gives  more  pleasing  encoded  images  to  the  human  viewer. 

Let  us  give  an  overview  of  this  coding  scheme  and  then  describe  the 
various  operations  in  detail.  The  block  diagram  of  this  scheme  is  given 
in  Figure  5.5.  First,  li  is  transformed  by  a  two-dimensional  unitary 
transform.  There  is  a  tradeoff  between  a  large  image  block  size  to 
extract  sufficient  edge  in format  ion  and  a  small  size  to  save  memory  and 
computat ional  effort.  We  find  that  a  good  compromise  is  achieved  by 
64  x  64  blocks.  Next,  applying  operator  on  V  (see  Fig.  5.5)  where  V 
is  the  transformed  image.  The  output  of  this  operation,  e,  is  segmented 
next  according  to  (5. 2. 1-3)  and  (5. 2. 1-4).  If  |c(itJ)|  >  t,  the  location 
(i.j)  is  recorded  as  one  of  the  edge  addresses,  otherwise  it  is  rejected. 


The  value  of  e  at  the  detected  edge  point  is  quantized  before  transmission 

to  the  receiver.  This  quantized  version  of  C,  e*,  is  also  needed  to 

calculate  V*  as  JT  *C*,  where  V*  is  viewed  as  the  deterministic  part  of 

the  transformed  image.  Subtracting  V*  from  V,  we  have  V  ,  the  stochastic 

a  8 

part  of  the  transformed  image.  A  conventional  transform  coding  scheme  is 

then  implemented  on  V^.  Finally,  V*,  c*  together  with  edge  address  map, 

arc  used  for  coding  and  transmission  over  the  communication  channel. 

At  the  receiver,  the  process  proceeds  in  a  fairly  obvious  way. 

First,  we  reconstruct  V*  based  on  the  received  e*  and  «5.  Then  summing 

d  ad 

It  with  V*,  yields  V*.  Finally,  the  inverse  unitary  transform  is  operated 
on  V*  to  recover  the  input  image. 
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We  now  describe  details  of  some  of  the  operations  that  need  to  be 
performed. 

1)  Run  Length  Coding 

The  purpose  of  this  operation  is  to  efficiently  transmit  the  edge 
address  map.  The  entropy  of  the  run  length  serves  as  a  measure  of 
bit  rate  required  to  recover  the  original  edge  address  map  without 
distortion.  For  higher  threshold,  t,  fewer  pixels  are  detected  as 
edRe  points.  In  this  situation,  in  order  to  prevent  excessively 
long  runs,  pseudo  runs  are  inserted  if  run  length  exceeds  a  fixed 
number  (64  in  our  experiment). 

2)  Fixed  Level  Quantization 

F.ssentially  we  have  only  one  quantizer  at  this  stage  and  e.(i,j)  is 

d 

assumed  to  have  a  I.aplacian  distribution  with  mean  zero  and  variance 

2 

a  which  has  to  be  calculated  for  each  value  of  t,  since  the  values 
Ld 

of  would  vary  with  t.  Due  to  the  fact  that  the  values  of 

C^(i,j)  are  usually  small,  it  has  been  found  that  a  two  bit  quantizer 

gives  a  very  accurate  reproduction  of  E  (i,j).  We  have  studied  the 

d 

effect  of  assigning  different  bits  to  the  quantizer.  Notice  that 
if  the  number  of  detected  edge  points  increases,  then  assigning  too 
many  bits  to  the  quantizer  would  mean  that  we  are  not  able  to  allocate 
sufficient  bits  to  the  transform  domain  quantizers,  and  this  would 
cause  the  objectionable  low  spatial  frequency  quantization  errors  in 
the  encoded  image. 

3)  Transform  Domain  Quantization 

| 

In  order  to  accomplish  the  task  of  data  compression,  we  must  allocate 
unequal  nuni>er  of  bits  to  the  quantizers  for  each  element  of  Vft.  The 
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fundamental  idea  of  transform  coding  is  that  high  energy  (activity) 
samples  are  quantized  more  accurately  than  low  energy  samples.  The 


variances  of  each  element  of  V  ,  the  stochastic  part  of  the  trans- 

s 

formed  image,  can  be  estimated  from  the  SDF  corresponding  to  the 
transform  used.  The  quantization  scheme  used  here  is  the  Lloyd-Max 
quantizer  with  the  assumption  that  each  element  of  Vr  is  Caussian 
in  distribution. 


From  the  above  equation,  clearly,  we  have  defined oT  as 


(5.3-2) 


Comparing  (5.3-2)  with  (5. 2. 1-4),  it  can  be  seen  that  X'  is  slightly 
difierent  f rom  tf,  the  way  presented  in  Fig.  5.5  is  easily  realizable. 

5)  Considerations  in  Deleting  Isolated  Point 

Edge  detection  is  often  affected  by  noise  contained  in  the  original 
image.  I f  we  examine  the  edge  picture  (picture  composed  of  the  detected 
edge  points  only)  closely,  see  Fig.  5.1,  we  would  find  many  isolated 
spots  spread  over  the  scene.  Actually,  we  have  no  sharp  distinction 
between  a  "true  edge  point"  or  a  "noise"  point.  A  simple  edge  detection 
scheme  in  the  presence  of  small  amount  of  noise  is  to  first  implement  the  edge 
extraction  algorithm,  then  followed  by  an  algorithm  that  removes  the  iso¬ 
lated  edges  (noise)  from  the  scene.  In  this  way,  one  could  get  a  much 


clearer  edge  profile.  We  have  Included  this  scheme  inside  edge 
detection  logic  box  of  Figure  5.5. 

An  advantage  of  this  transform  coding  method  is  that  it  gives  us  a 
flexible  way  to  meet  any  requirement  of  getting  a  good  visual  image  sub¬ 
ject  to  a  fixed  bit  rate.  There  are  three  adjustable  parameters  associated 
with  the  entire  process,  viz, 

a-  F.dge  point  selection  threshold 

The  basic  formula  is  given  by  t  *  CO,  where  0  is  the  rms 
value  of  C ( 1 , j ) ,  and  c  is  an  adjustable  constant.  When  c  “  no 
pixel  will  be  detected  as  an  edge  point.  Some  of  the  operations 
in  Figure  5.5  could  be  bypassed  in  this  case  and  a  conventional 
transform  coding  method  is  directly  performed  on  the  image  0. 

b.  Edge  quantizer’s  bit  (n^  bits/edge) 

We  have  made  experiments  concerning  this  fixed  integer  hit's 
effect  on  a  bit  allocation  strategy.  The  various  bit  rates  used 
are  1,  2,  3  and  4.  The  results  will  be  presented  in  the  next 
section.  Experimentally,  it  has  been  found  that  2  bits  (4  levels) 
gives  a  proper  bit  selection. 

c.  Desired  bit  rate  (n) 

It  is  the  desired  overall  average  bit  rate  in  bits/pixel  in 
the  final  encoded  image. 

The  above  mentioned  parameters  are  related  by 

n  -  n^  +  n£  +  n^  (5.3-3) 
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n ^  :  Average  bit  rate  for  transform  coding  of 


n 


2 


n3 


Average  bit  rate  for  coding  the  amplitudes  of  detected 
edge  points 

n  x  (Number  of  Detected  Edge  Points) 

— - -  (5.3-4) 

(Number  of  Image  Pixels) 

Bit  rate  required  to  recover  the  edge  address  map,  jL. 

(Entropy  of  *9  )  x  (Number  of  Detected  Edge  Points) 

. — - -  - -  (5.3-5) 

(Number  of  Image  Pixels) 


Once  the  edge  point  selection  threshold  is  given,  the  number  of 

detected  edge  points  is  automatically  determined.  For  a  fixed  n  ,  n„  and 

e  2 

n}  can  be  calculated  via  (5.3-4)  and  (5.3-5),  respectively.  I f  we  desire 

g 

to  have  bit  rate  of  n  (or  -  data  compression)  in  the  final  encoded  image, 
n^  has  to  be  calculated  via  (5.3-3). 


5.4  Experimental  Results  and  Comparisons 

A  series  of  computer  simulations  has  been  conducted  to  evaluate  the 
performance  of  this  scheme  called  Feature  transform  coding.  A  256  x  256 
image  was  first  divided  into  16  equal-sized  image  blocks,  each  of  size 
64  *  64.  Each  image  block  was  feature  transform  coded  independently.  In 
all  the  experiments.  Cosine  transform  based  NCI  model  was  used  and  several 
definitions  of  performance  measurements,  i.e.,  SNR,  N>iSE,  etc.,  are  defined 
in  section  3.5  of  Chapter  3.  Two  inages  of  different  category  (statistics) 
were  used  in  order  to  test  how  sensitive  the  scheme  was  to  the  image 
content . 

The  quantities  np  and  t  are  the  major  parameters  in  the  design  of 
Feature  transform  coding.  In  our  study,  it  is  done  experimentally  as 
presented  in  the  next  section. 
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5.4.1  Feature  Transform  Coding 

In  order  to  find  the  proper  choice  of  n  (edge  quantizer's  hit)  for 

e 

different  values  of  t  (threshol d) ,  we  plot  the  curves  of  SNR  and  NMSE 

versus  various  settings  of  n^  for  t  “  1.5a  and  t  ”  2o.  Figures  5.6-5. 7 

and  Figures  5. 10-5. 11  are  the  results  for  Girl  and  Chemical  Plant  images, 

respectively.  It  can  he  seen  from  these  figures  that  n  must  be  equal  to 

e 

2,  since  it  gives  the  minimum  NMSF.  or  maximum  SNR  of  all  the  bit  rates. 
After  determining  n  ,  we  also  plot  the  curves  of  SNR  and  NMSE  versus 
various  thresholds.  Figure  5.8  demonstrates  the  curve  for  the  Girl  image, 
and  Figure  5.12  is  for  the  Chemical  Plant  image.  It  is  observed  that  the 
SNR  increases,  while  the  threshold  decreases  for  moderate  bit  rate  (1~1.5) 
of  both  images.  For  high  bit  rate  (*2)  this  phenomenon  is  still  observed 
for  Chemical  Plant  image,  but  not  for  Girl  image.  This  could  be  because 
the  Chemical  Plant  image  is  a  complicated  scene  in  which  more  edge  points 
must  Be  separated  from  the  background  and  encoded  adequately.  For  low 
bit  rate  («0.5),  as  expected,  this  scheme  performs  slightly  better  than 
the  standard  transform  coding  (in  our  case  NCI  Cosine  transform  coding) 
for  higher  thresholds  only.  Based  on  these  observations,  it  seems  that 
t  ■  1 . 5<J  and  n^  ■  2  serve  as  the  best  parameters  to  he  used  to  encode  any 
given  image  at  rate  greater  than  0.5  bit  /pixel. 

The  comparisons  of  Feature  transform  encoded  and  standard  transform 
encoded  images  are  shown  in  Figure  5.14  for  the  Girl  image  and  correspond¬ 
ing  ones  for  the  Chemical  Plant  image  are  shown  in  Figure  5.15.  In  all 
the  images  the  Feature  transform  coder  produces  a  "better"  appearing  image 
(at  bit  rate  >  0.5)  than  the  standard  transform  coder.  It  is  easily  seen 
that  there  is  a  significant  difference  around  the  Girl's  face  and  at  the 
sides  of  the  freeway  in  the  Chemical  Plant  image.  In  general,  the  edges 
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are  more  clearly  defined  on  the  Feature  transform  encoded  image.  The 
above  discussion  is  also  supported  by  referring  to  the  coding  results 
shown  in  Tables  5.2  through  5.9,  which  give  quantitative  comparisons  of 
the  above  two  mentioned  schemes.  For  easy  reference,  the  performances  of 
the  standard  transform  coder  at  different  bit  rates  are  listed  along  with 
the  Feature  transform  coder.  The  p -rformance  improves  by  about  1.5  dB 
at  1  bit  rate  for  both  images  and  by  2dB  at  2  bit  rate  for  Chemical  Plant 
only.  It  is  of  interest  to  note  that  at  2  bit  rate  Figure  5.1A(f)  has 
less  SNR  value  than  Figure  5.14(c).  However,  by  directly  looking  at  the 
pictures,  one  can  say  that  Figure  5.14(f)  gives  a  Letter  appearing  image 
since  the  edges  are  more  sharp  regardless  of  the  contrast  of  the  background. 

5.4.2  Classification  Feature  Transform  Coding 

In  this  method  the  edge  inform.it  ion  map  <9  and  edge  point  r.  ,  are 

a  d 

encoded  and  transmitted  as  before,  except  that  the  stochastic  portion  of 
the  transformed  image,  V  ,  is  classified  as  belonging  to  one  of  K  pre¬ 
determined  classes.  In  our  case,  K  ■  2  and  each  class  ip.  assumed  to  be 
equiprobab! e.  The  class i f i cat  ion  maps  for  the  Girl  and  Chemical  Plant 
images  are  shown  in  Fig.  4.11(a)  and  Fig.  4.11(b).  The  block  diagram  of 
this  scheme  is  analogous  to  Figure  5.5.  Howe\?r,  the  Transform  Domain 
Quantization  box  now  contains  two  types  of  hit  allocation  maps,  one  for 
each  class. 

This  method  gives  an  additional  improvement  of  SNR  over  the  Feature 
transform  coder  about  1  dB  at  hit  rate  >  1.5  for  the  Girl  image.  But  it 
failed  for  the  Chemical  Plant  image.  Figures  5.9  and  5.13  show  the  coding 
r-  ..Its  versus  various  thresholds  with  np  *  2  assigned  to  the  edge  point 
•  >  it  tzar.  The  reason  for  the  discrepancy  could  be  that  the  Chemical  Plant 


159 


Image  is  a  complicated  scene,  and  once  enough  detected  edge  point  infor¬ 
mation  has  been  subtracted  from  the  original  nonstationary  image  field,  the 
resultant  image  field,  i.e.  Vg,  would  tend  to  have  equal  variances  for 
each  image  block,  i.e.,  the  corresponding  pixels  of  the  remaining  field 
Vg  of  each  image  block  have  approximately  equal  variances.  However,  this 
kind  of  property  is  violated  for  a  simple  structure  image,  i.e..  Girl 
image,  where  after  the  edge  point  information  is  subtracted,  some  refined 
work,  for  example  classification,  could  be  applied  on  if  the  ultimate 
goal  of  coding  scheme  is  to  get  higher  SNR  in  the  final  encoded  image. 

We  have  shown  the  encoded  images  of  this  scheme  along  with  Feature  trans¬ 
form  encoded  Images  In  Figures  5.16  and  5.17. 

5.4.3  Conclusions 

Based  on  the  above  experiments  and  the  theory  discussed  in  Chapter  4, 
the  following  conclusions  are  made: 

(1)  Regardless  of  the  computational  complexity.  Feature  transform  coding 
scheme  can  improve  the  signal  to  noise  ratio  over  the  standard  trans¬ 
form  coding  scheme  (based  on  NCI  Co9ine  transform  model)  by  1  to 

2  dB  at  hit  rates  above  1  bit/pixel  for  usual  images. 

(2)  The  complexity  of  Feature  transform  coding  scheme  Is  simpler  than  the 

technique  proposed  by  Yen  et  al  [  74  1  where  they  ignored  the  problem 

of  computation  tine,  complexity,  and  storage.  Meanwhile,  they  had 

more  linage  dependent  parameters  that  need  to  he  identified.  In  our 

case,  onlv  two  parameters,  i.e.,  n  and  t,  are  sufficient. 

e 

(3)  The  Feature  transform  coding  scheme  is  useful  when  it  is  desired  to 
maintain  the  attractive  characteristics  of  the  transform  coding 
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scheme,  i.e.,  less  sensitivity  to  the  image  statistics  (scene  to 
scene  variation),  as  well  as  to  incorporate  the  physiological  reac¬ 
tion  of  the  human  visual  system  in  the  overall  communication  system 
design. 

(4)  The  Classification  Feature  transform  coding  scheme  is  suited  for 

simple  structured  images  at  bit  rate  >  1  bit/pixel.  The  complexity 
over  the  Feature  transform  code  is  marginal  and  its  performance  is 
improved  by  reducing  the  sire  of  the  image  block  or,  in  other  words, 
the  number  of  image  blocks  is  increased. 


t 
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(Edge  Addresses) 


Figure  5.2  Realization  of  Stochastic  Decomposition 


(c)  Original  (d)  Edge  Pixels 


Figure  5.31  Image  Decomposition  of  the  Girl  Image  with  3111  (t  »  2o) 
Detected  Edge  Points 
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(c)  Original 


(d)  Edge  Pixels 


Figure  5.41  laage  Decomposition  of  the  Chemical  Plant  Image  with  3653 
(t  •  2a)  Detected  Edge  Points 
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(e)  U  ♦  2.5U.  U  ♦  JU., 

s  d  8  d 

Figure  5 . 4 T I  Image  Enhancement  of  the  Chemical  Plant  Image  with  3653 

(t  ■  2 n)  Detected  F.dge  Points 
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y  (Edge  Addresses) 


Figure  5.5  Transform  Coding  Based  on  Stochastic  Segmentation 
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Figure  5.7  Feature  Transform  Coding  Result.  for  the  Girl  Image  ver.u.  Varlou*  Edge  Quantizer's 
Sit.  with  Jill  (t *2o)  Detected  Edge  Point. 
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Figure  58  Feature  Transform  Coding  Result*  for  the  Cirl  Image  versus  Various  Thresholds 
with  ?  Bits  Assigned  to  the  Edge  Quantiser 
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Figure  5.S  Classification  feature  Transform  Codine,  Results  for  the  Cl  rl '  Image  versus 
Various  Threshold*  with  7  Bit*  Assigned  to  the  Edge  Quantiser 


170 


Signoi  to  Noiso  Rc*o  (dB)  2  Signo!  to  Noise  Retie  (dB) 


32 


31 


30 


29 


26 


Edge  Ouontirer’s  Bit 


Tronvform 

Coding 


gure  5.10  Feature  Transform  Coding  itosulta  for  the  Chemical  Plant  Image  vrrtui  Various  Edge 
Quentixrr'*  Bit*  with  7079  (t  -  1.5a)  Detected  Edge  Point* 
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Figure  5.11  Feature  Transform  Cottier.  Be*tilt«  for  the  Chemical  Mam  Image  ver»u*  Various  Edge 
Quantiser's  Bits  with  J651  (t  "  ?")  Detected  Edf.e  Point* 
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(a)  1  Bit /Pixel , 

Standard  Trans  form 


(d)  1  Bit/Pixel, 

Feature  Transform 


(e)  1.5  Bits/Pixel, 
Feature  Transform 


(b)  1.5  Bits/Pixel, 

Standard  Transform 


(c)  2  Blta/Plxel,  (f)  2  Bits/Pixel, 

Standard  Transform  Feature  Transform 


Flgurp  5.14  Comparison  of  Standard  Transform  F.ncoded  and  Feature 
Transform  F.ncoded  (n  •  2  Bits,  t«1.5o,  5737  Detected 
EdRe  Points)  Girl  ImSftes 
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(d)  1  Bit/Pixel, 

Feature  Transform 


(a)  1  Bit/Pixel 

Standard  Transform 


(b)  1.5  Bits/Pixel, 

Standard  Transform 


(e)  1.5  Bits/Pixel 

Feature  Transform 


(c)  2  Bits/Pixel ,  (f)  2  Bits/Pixel, 

Standard  Transform  Feature  Transform 

Figure  5.15  Comparison  of  Standard  Transform  Encoded  and  Feature  Transform 
Encoded  (n  •  2  Bits,  t«1.5o,  7079  Detected  Edge  Points) 
Chemical  Plant  Images 
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(a)  l  Bit/Pixel 


(d)  1  Bit/Pixel, 
Classification 


(e)  1.5  Blts/Plxel 
Classification 


(c)  2  Bits/Pixel 


(f)  2  Blts/Plxel, 
Classification 


Figure  5.16  Feature  Transform  Encoded  Girl  Image  with  n 
f  •  1.5o  (5717  Detected  Edge  Points) 
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Figure  5.17  Feature  Transform  Encoded  Chemical  Plant  Image  with  n 
and  t  •  I.Vj  (7079  Detected  Edge  Points) 
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CHAPTER  SIX 


SUMMARY  AND  RECOMMENDATIONS  FOR  FUTURE  RESEARCH 

This  thesis  has  shown  that  the  application  of  the  stochastic 
laage  models  based  on  PDEs  to  the  problem  of  image  data  compression 
has  achieved  considerable  success.  The  use  of  image  models  has 
shown  the  connection  between  predictive,  hybrid  and  transform  coding 
schemes.  From  the  choice  of  an  approximate  Image  model  which  gives  the 
best  approximation  to  the  actual  statistics  of  a  class  of  images 
being  processed,  on  would  be  able  to  determine  the  most  efficient 
data  compression  scheme. 

Simple  basis  restriction  error  experiments  were  performed  to 
compare  the  energy  compaction  property  of  various  image  models.  The 
results  show  that  both  the  noncausal  and  semlrausal  models  give 
superior  performance  over  the  causal  models.  F.vidently,  different 
spectral  density  function  shapes  of  image  models  were  the  dominating 
factor  in  this  aspect  of  study. 

The  performance  of  the  hybrid  coding  scheme  was  significantly 
improved  by  two  very  simple  methods,  i.e.,  (1)  Adaptive  Variance 
Estimation:  Adapting  the  variance  of  the  prediction  error,  in  each 
DPCM  channel,  the  updated  variance  was  used  to  adjust  the  spacing  of 
the  quantizer  levels.  (2)  Adaptive  Classification:  Adapting  the  bit 
rate  to  the  local  variance  of  each  Image  column. 

The  adaptive  variance  estimation  scheme  was  shown  having  the 
potential  of  maintaining  the  performance  of  a  hybrid  coder,  designed 
for  a  nominal  statistics,  in  the  face  of  changing  statistics.  The 
adaptive  classification  scheme  performed  more  effectively  at  low  bit 
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rate  and  was  robust  with  respect  to  channel  errors.  Hybrid  coding  of 
noisy  images  was  achieved  by  replacing  the  predictor  in  the  DPCM  loop 
by  a  Kalman  filter.  Experimental  results  showed  that  the  two  bit  rate 
quantizer  gave  a  restored  image  with  SNR  slightly  less  than  the 
infinite  rate  (without  quantizer  in  the  DPCM  loop). 

A  new  technique  called  feature  transform  coding  for  efficient 
transform  coding  of  images  was  presented.  The  experimental  results 
of  applying  this  technique  to  actual  image  data  (Cirl  and  Chemical 
Plant  images)  show  that  it  provides  us  a  flexible  way  of  getting  a 
good  visual  image  (subject  to  a  fixed  bit  rate)  while  maintaining  the 
mean  square  error  performance.  It  is  evidenced  by  referring  to  Figures 
5.6  -  5.13.  The  incorporation  of  classification,  a  scheme  analogous 
to  the  adaptive  classification  hybrid  coding,  results  in  further 
improvement  in  the  performance  for  simple  structured  images  at 
bit  rate  1  bit/pixel.  Moreover,  its  complexity  over  the  feature 
transform  coder  is  marginal. 

From  a  practical  point  of  view,  the  most  important  results  of 
this  research  is  that  it  provides  a  theoretical  framework  for 
optimizing  the  design  of  an  intraframe  coding  system  for  two  dimensional 
data  modeled  by  random  fields.  The  various  coding  algorithms  developed 
here  can  be  easily  implemented  for  designing  practical  coding  systems 
for  images. 

Areas  in  need  of  further  investigation  are: 
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1)  Inclusion  of  NC2  and  NC3  models  to  the  design  of  feature 
transform  coder.  Since  only  NCI  model  based  feature  transform 
coder  has  been  fully  explored. 

2)  Extension  to  the  higher  order  stochastic  image  models  which 
might  be  needed  to  represent  image  spectral  density  functions 
(SDKs)  more  accurately. 

3)  Application  of  data  compression  schemes  to  biomedical  images. 

4)  Study  the  effects  of  Incorporating  variance  estimation  scheme 
in  the  feature  transform  coder. 

5)  Use  of  stochastic  models  in  the  design  of  coders  in  the  presence 
of  channel  errors. 
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APPENDIX  A 

QUANTIZER  CONSIDERATIONS 

A.l  Introduction 

The  basic  principle  in  designing  a  quantizer  is  to  determine  the 
optimum  quantization  levels  so  as  to  minimize  the  overall  mean  square 
error  between  the  original  signal  and  quantized  signal.  The  exact 
approach  to  quantization  problem  can  only  be  solved  by  wasting  a  con¬ 
siderable  amount  of  computing  time  and  gives  little  insight  for  someone 
who  has  a  new  quantization  problem  and  needs  it  to  be  solved  in  a  small 
period  of  time.  The  useful  approximation  to  the  optimum  quantizer  was 
first  proposed  by  Algazi  |  2  ].  In  this  Appendix,  we  have  modified  his 
formulas  slightly  and  tabulated  the  characteristics  of  the  resulting 
quantizer  (so  called  compandor)  for  Gaussian  and  l.aplaeian  distributed 
signals  together  with  the  optimum  nonuniform  (in  the  minimum  mean  square 
error  sense)  and  optimum  uniform  quantizers  which  have  been  considered 
by  Max  |  44  )  and  others.  Also  rate-distortion  functions  f(n)  are  listed 
which  were  obtained  by  a  piecewise  exponential  curve  fitting  technique  to 
the  tabulated  quantizer  numerical  results. 

A. 2  Quantizer  Structures 

A  typical  quantizer  characteristic  is  shown  in  Fig.  A.l.  The  quan¬ 
tizer  output  y  is  assigned  to  one  of  N  discrete  values  (Fig.  A. 1(a)  and 
Fig.  A.l(h)  are  for  N  even  and  odd  respectively)  depending  on  the  value 
of  input  x.  The  x^  and  are  called  the  "decision  levels"  and  "quantlza- 
X  tion  levels"  respectively.  If  the  probability  density  function,  p(x),  of 

the  quantizer  input  is  given  and  is  assumed  to  be  symmetric,  one  would 
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have  different  quantizers  according  to  the  diverse  assumptions  of  design 
constraints.  In  this  section,  we  briefly  review  some  of  the  useful 
quantizer  reults.  The  interested  reader  can  find  detailed  information  in 
references  12, 40, **4, 51, 72] 

A. 2.1  Optimum  Nonuniform  Quantizer 

In  this  quantizer,  the  design  procedure  is  to  find  the  N  y^'s  and 

associated  x^'s  in  Fig.  A.l  so  as  to  minimize  the  distortion,  D,  which  is 

defined  as  the  expected  value  of  W(e),  where  W  is  an  error  weighting 

function  and  e  is  the  quantization  error.  For  simplicity,  W(r)  is  assumed 
2 

to  be  c  ,  the  square  value  of  the  quantization  error.  The  distortion  is 
obtained  as 

Case  1:  N  •  even 


where  y  ■  0,  Xj  >  0,  xN+1  -  ®,  and  lies  between  x^j  and  . 

y- 

In  order  to  minimize  D  for  a  fixed  N,  the  necessary  conditions  are 
obtained  by  different  iating  Dwith  respect  to  the  x^'s  and  y^  's  and  setting 
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the  results  equal  to  zero.  After  simple  numerical  manipulation,  we  have 


Case  1:  N  »  even 


Xj  “  1  “  (A. 2.1- 3a) 

X1 

(yj~x)p (x)dx  -  0  1  -  l,...,y  (A.2.1-3b) 


Case  2:  N  ■  odd 

xA  -  *sCy t 1+1 )  *  •  l,...r^r  (A.2.1-4a) 

(yi-x)p(x)dx  -  0  1  ■  2,...,^—  (A.2.1-4b) 

The  above  equations  moan  that  x^  lies  midway  between  the  adjacent  quanti¬ 
zation  levels  and  y^  Is  the  centroid  of  the  area  under  p(x)  and  lies 

between  the  adjacent  decision  levels. 

A. 2. 2  Optimum  Uniform  Quantizer 

This  quantizer  Rives  the  minimum  distortion  under  the  constraint 

"  yj-y^i  ■  q  v  l  (a. 2.2-1) 

where  q  is  the  equal  length  spacing  between  adjacent  decision  levels  or 
quantization  levels. 

This  is  an  easier  problem  to  be  solved  than  the  previous  optimum  non- 
uniform  quantizer.  For  this  case,  the  distortion,  D,  becomes  a  function 
of  q  and  of  any  particular  value  of  N.  If  the  probability  density  function 
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of  the  quantizer  Input  is  assumed  to  be  known  and  symmetric,  then  a 
symmetric  answer  for  the  decision  levels  and  quantization  levels  will  be 
obtained.  The  distortion,  D,  can  be  obtained  as 


Case  1:  N  -  even 


N  1  4 

t  -1  iq 


2  l  |  [(^2^)q-*)  p<x)dx  * 2  0r^)q~x]  p(x>dx 

1"l(1-1)q  (!-l)q 


(A. 2. 2-2) 


Case  2:  N  ■  odd 


(i+*s)q 


f  2  t  N  f  r  T 

D  -  |  x2p(x)dx  +  2  y  |  (iq-x) 2p(x)dx  +  2j  l(^2^)q_XJ  P^dx 

-f  '"a-vi  (£*)<.  (*. 


(A. 2. 2-3) 


The  optimum  uniform  quantizer  step  size  q  can  be  obtained  by  differentiat¬ 
ing  D  with  respect  to  q  and  setting  the  results  equal  to  zero.  We  have 


Case  1 :  N  -  even 


-  1  iq  » 

®  .  l\  (21-1)  f  [(^-)q-^  P(.W.  ♦2C»-1>[  gS|L)q-*]p<«><*  - 

«-»q  (».,), 


(A. 2. 2-4) 


Case  2:  N  “  odd 


(  i+*j)q 


2  '•■•'I  •>  « 

^  -  2  l  i  j  (iq-x)p(x)dx  +  (N-l)  QT^)q"xJ  p*x)dx  -  0 

1-1  (i-S)q  (Hl2)_ 

i  2 


(A. 2. 2-5) 
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In  either  case,  the  problem  is  quite  acceptable  to  machine  compu¬ 


tation  when  p(x)  and  N  are  specified. 

A. 2. 3  Useful  Approximation  to  the  Optimum  Nonunlform  Quantizer  (  2  ] 

Since  a  quantizer  is  often  easier  to  implement,  we  consider  the 
nonunlform  quantizer  hh  the  case  cade  of  two  nonlinear  devices  and  of 
a  uniform  quantizer  as  shown  in  Fig.  A. 2.  With  a  given  uniform  quantizer 
and  a  given  probability  density  function,  p(x),  the  tvo  nonlinear  devices 
f(-)  and  g(-)  in  Fig.  A. 2  are  chosen  so  as  to  minimize  the  distortion, 

D  •  E[W{ g(y)-x) ] ,  where  W  as  defined  in  section  A. 2.1  is  an  error  weight¬ 
ing  function.  The  distortion  is  made  up  of  two  parts: 

1.  The  distortion  obtained  within  the  range  of  the  uniform 
quantizer. 

2.  The  distortion  in  the  tails  D^.  for  which  the  nonlinear  devices 
f(*)  and  g(')  are  completely  ineffective. 

The  analog  system  is  shown  in  Fig.  A. 3.  The  quantization  noise  due 
to  uniform  quantizer  is  assumed  to  be  independent  of  the  input  signal  for 
a  large  number  of  quantization  levels,  the  effect  of  which  has  been  dis¬ 
cussed  in  detail  by  many  researchers.  Referring  to  Fig.  A. 3,  the  two 
extreme  points  of  the  output  of  the  first  nonlinear  device  is  a  dominating 
factor  in  the  determination  of  the  overall  distortion,  in  other  words,  D 
is  a  function  of  f(x  )  and  f(x  .  ),  and  f(x  )-f(x  .  )  «  Nq,  q  being 
the  step  size  of  the  tmiforra  quantizer.  The  output  of  the  second  non¬ 
linear  device,  g(y),  is  expressed  in  terns  of  n  and  f(x)  as 

R<y)  •  ♦  n]  -  g|f(x)]  +  ng'If(x)] 
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(A.2.3-1) 


If  Che  quantisation  error  Is  snail,  a  Taylor  series  expression  of 


g(y)  gives  the  above  equation  directly.  And  the  expression  for  Dn  can 
be  written  as 

3.  x 
f 2  t  max 

D  -  W(q(f(a)+6]  -  a}p  (a)p  (S)dadB  (A. 2. 3-2) 

n  J  q  )  x  n 

"  2  Xmin 

The  limits  on  the  Integrals  are  — ^  to  ^  for  the  integral  in  8  and 

x  to  x  (f(x  )  _  *(*...(„>  *  N<l)  for  the  integral  in  a. 
tain  max  max  ram 

2 

If  we  assume  W (c)  -  e  and  g(y)  is  sufficiently  smooth  and  the 
quantization  error  is  small,  then  following  the  development  in  [  2  ],  one 
would  be  able  to  obtain 


Dn  ’  , 

12N2 


max  1 

(px(a)  ]Tda 


min 


1  3 


(A. 2. 3-3) 


If  we  take  into  account  the  probability  of  occurence  of  and  I)T, 
we  have 


D-D 


—  X  ®  -i  X 

f  min  r  I  r  min 

1  -  I  px(a)dn  -  I  px(a)dn  +  j  [ 

-  x  J  -*> 


max 


f 80 

♦  [(X  -|)]2p,(a)do 

J  max  *  * 


((x»in  +  !>  ‘  xl  px(a)d“ 


(A. 2. 3-4) 


For  a  symmetric  probability  density  function,  we  set  x^^  “  “\jn  *  3 


and 


3N 


3r 


j  [px(a)Pda  1— 2 J  px(a)dn  +  2  j  ( (X-|)-x]2px<a)<ta  (A. 2. 3-5) 


Equation  (A. 2. 3-5)  has  to  be  minimized  by  proper  choice  of  X.  Once 
X  is  determined ,  then  it  can  be  shown  that 

x 

Nqj  [px(a)Fda 

f(x)  -  — - - .  (A. 2. 3-6) 

2j  [px(a)]7da 
0 

It  i9  noted  from  the  above  equation  that  the  step  size  of  the  uniform 
quantizer  will  not  affect  the  resulting  nonuniform  quantizer. 

The  characteristic  of  this  quantizer  for  2m  (m  ■  1,...,9)  quantiza¬ 
tion  levels  is  given  and  discussed  in  the  next  section.  We  mention  here 
that  the  derivation  given  in  (  2  ]  erred  in  assuming  the  two  extremes  of 
the  uniform  quantizer  shown  in  Fig.  A.  3  to  be  ±  (^2^) ^  *n8tea<*  4^Jq* 

A. 3  Results  and  Comparisons 

The  structures  for  the  optimum  uniform,  optimum  nonuniform  and 
approximate  nonuniform  (compandor)  quantizers  are  described  in  the  pre¬ 
vious  sections.  In  this  section,  we  present  some  nice  characteristics  of 
the  approximate  nonunifora  quantizer  along  with  the  optimum  uniform  and 
optimum  nonuni  form  quantizers.  The  equations  defining  the  optimum  non- 
uniform  quantizer  yield  closed-form  solution  only  in  certain  special  cases, 
depending  solely  on  the  probability  density  function  of  the  quantizer  in¬ 
put  signal.  For  example,  even  for  the  Caussian  distributed  signal,  for 
lack  of  the  closed-form  solution  for  the  x^  (decision  levels)  and  y^  (re¬ 
construction  levels),  some  sort  of  numerical  solution  must  be  sought. 

Max  [  44  ]  suggested  an  iterative  procedure  to  calculate  x^  and  y^  and 
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tabulated  the  values  of  x^  and  for  N  (number  of  levels)  up  to  36,  since 
the  time  required  to  determine  the  quantizer  structure  increases  exponen¬ 
tially  with  the  total  number  of  levels  of  the  quantizer.  Max's  [44  1  work 
cannot  be  extended  beyond  N  ■  36.  Kurtenbach  and  Wintz  (40  )  applied  the 
same  numerical  approach  and  tabulated  the  value  of  xJ  and  yJ  to  N  «  512 
(N  “  2n,  n  •  9,  n  :  number  of  quantizer  bits)  for  the  Gaussian  distributed 
signal.  They  confessed  also  that  the  determination  of  the  quantiser 
structure  for  m  ■  6,  7,  8,  9  is  limited  because  of  the  computer  time  re¬ 
quired  for  each  case.  Paez  and  Glisson  (51  ]  published  the  quantizer 
results  for  the  Laplacian  and  Gamma  distributed  signals  facing  the  same 
kind  of  difficulty,  and  tabulated  the  data  only  up  to  N  ■  32.  So  the 
exact  approach  to  the  optimum  nonuniform  quantizer  gives  little  insight 
and  sometimes  prsents  someone  who  has  a  new  quantization  problem  with  a 
considerable  amount  of  computation  time  to  obtain  a  significant  result. 

The  second  advantage  of  designing  the  approximate  nonuniform  quantizer  is 
that  it  is  easily  iroplementable  with  only  a  marginal  increase  in  its 
complexity  compared  to  the  optimum  uniform  quantizer. 

Tables  A.l  and  A. 2  list  the  approximate  average  amount  of  CDC  CYBER 
173  computer  time  required  to  determine  the  quantizer  structure  for  the 
optimum  nonuniform,  optimum  uniform  and  approximate  nonuniform  quantizers. 

A  comparison  of  these  three  quantizers  can  be  easily  made  from  the  entries 
of  Tables  A.l  and  A. 2  to  conclude  that  the  approximate  nonuniform  quantizer 
is  the  best  choice  as  far  as  computation  time  is  concerned. 

Note  that  we  have  used  our  own  numerical  approach  to  solve  the 
values  of  x^  and  y^  for  a  given  quantizer.  The  procedure  is  described 
as  follows: 
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If  the  probability  density  function  is  known,  we  make  an  initial 
guess  for  Xj  (I  ■  l,...,^  -  1),  where  the  x^  are  uniformly  spaced. 

2 

is  assumed  to  be  equal  to  five  times  the  standard  deviation  of  the  input 
signal.  Then  yA  is  calculated  (i  -  l,...,y)  by  solving  (A.2.1-3b).  x^^ 

(i  ■  1 , . . .  -  1)  is  updated,  denoted  as  x^,  by  using  Newton's  numerical 

method  on  (A.2.1-3a)  compute  resulting  error  e  •  -  x^.  If  it  satisfies 

the  required  error  bound,  then  x^  was  chosen  correctly,  otherwise  the 
process  is  repeated.  Fig.  A. 4  shows  the  flow  chart  of  this  numerical 
procedure . 

The  correctness  of  the  computer  program  was  checked  by  printing  out 
the  numerical  value  of  each  data  entry  to  five  digits  after  the  decimal 
point.  Comparing  the  results  with  Max's  [  44  ] ,  we  have  precise  value 
for  each  x^  and  y^,  and  furthermore  our  accuracy  is  extended  to  five 
digits  instead  of  four  digits  after  the  decimal  point. 

As  a  by-product  of  determining  the  quantizer  structures,  most  of  the 
useful  information  quantities  were  obtained.  There  are  (a)  the  mean  square 
error  (M.S.E.)  introduced  in  the  quantization  process,  (b)  the  entropy  at 
the  quantizer  output,  (c)  the  signal  to  noise  ratio,  (d)  the  rate  distoi>- 
tion  function  or  source  information  rate.  In  addition,  all  the  data 
entries  tabulated  are  constructed  for  quantizer  input  signal  with  unit 
standard  deviation  (r.m.a.  value).  To  obtain  the  'decision  levels'  x^ 
or  "quantization  levels'  y^  for  signals  with  standard  deviation  different 
from  unity,  we  simply  need  multiply  the  given  values  of  data  entries  by 
the  actual  standard  deviation.  The  M.S.E.  is  found  by  multiplying  the 
given  M.S.E.  by  the  variance  of  the  signal. 
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In  all  the  measurements,  after  obtaining  the  values  of  and  y^, 
the  mean  square  errors  (M.S.F..)  are  calculated  by  (A. 2. 1-1),  (A. 2. 2-2) 
and  (A. 2. 3-5)  for  optimum  nonuniform,  optimum  uniform  and  approximate 
nonuniform  quantizer  respectively.  Here  we  assume  that  the  distribution 
of  the  quantizer  input  signal  and  quantizer  structure  is  symmetric  and 
that  the  total  number  of  quantizer  output  levels  is  even.  This  constraint 
is  not  restrictive.  For  example,  we  have  plotted  Tables  A.10-A.15  for 
N  »  2,.. .,36,  which  has  both  cases  N  “  even  and  N  ■  odd. 

For  entropy  it  is  possible  to  treat  the  quantizer  output  as  a  dis¬ 
crete  data  source  and  can  be  calculated  from  the  following  formula. 


Entropy  ■  -2  ]  p.  log-  p  bits  per  sample 

i-1  1  1  1 


where  pt  -  p  ly-y ± 1  -  p  [x 


f*1 

i-i  -  *  i  '  J 


p(x)dx  i  *•  1 


(A. 3-1) 

N 

•2  * 


Gaussian  and  Laplaolan  signal  distributions  were  considered.  These 
are  the  two  commonly  used  assumptions  for  the  signal  distribution  in  the 
image  processing  area. 


Guassian  density: 


,  v  1  "2a* 

pM '  ms  e 


(A. 3-2) 


rms  value  ■  a 
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The  exact  formula  for  In  either  case  can  be  expressed  as  follows: 


Laplacian  distributed: 


>t  -  f‘  P(.)dx  -  f‘  |  «-ai'tldx  -  (o'*'1  -  1-1 . H 


vi-l 


i-1 


(A.  3-4) 


Gaussian  distributed: 

x. 


3( 

t  “  |  p(x)dx  -  |  — e  2a  dx-sTerf^//?©)  -  erf (x^/ZSa^  * 


i-1 


i-1 


i  ■  1 » • • • »2 


-  (A. 3-5) 


The  signal  to  noise  ratio  (S.N.R.)  is  defined  as 


S.N.R.(dB)  -  10  log,- - .  (A. 3-6) 

10  mean  square  error  (M.S.E.)  ' 


where  we  assume  unity  variance  for  the  quantizer  input  signal. 

Tables  A. 3  -  A. 8  tabulate  the  characteristics  of  three  different 
quantizers  with  Gaussian  and  l.aplacian  distributed  quantizer  input  signal. 

In  order  to  better  compare  the  performances  of  these  three  quantizers, 
we  have  plotted  the  curves  of  all  three  quantizers  into  one  single  graph 
according  to  differene  categories,  i.e.,  mean  square  error,  signal  to  noise 
ratio,  and  entropy.  These  are  shown  in  Figs.  A. 5  -  A. 10.  It  is  seen 


t  We  use  the  National  Bureau  of  Standards  definition  for  the  error 
function,  i.e.. 


erf (x)  -  JL 

o 
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from  these  figures  that  approximate  nonuni  form  quantizers  (compandor) 
perform  virtually  as  good  as  optimum  nonuni  form  quantizers,  since  the 
characteristic  curves  of  these  two  quantizers  are  so  close  together. 

Fig.  A. 11  shows  the  unequal  probabilities  of  quantizer  output  levels  for 
the  case  of  N  ■  8.  Again,  the  similarity  of  the  compandor  and  the  non- 
uniform  quantizer  can  be  easily  observed.  Another  significant  result  that 
can  be  observed  from  Figs.  A. 5  -  A. 10  is  that  the  advantage  of  the  optimum 
nonuniform  quantizer  and  approximate  nonuniform  quantizer  (compandor)  over 
the  optimum  uniform  quantizer  starts  with  n  ■  2  and  increases  as  n  in¬ 
creases.  This  effect  is  due  to  the  fact  that  both  former  quantizer  struc¬ 
tures  have  more  quantization  levels  near  the  origin  than  the  optimum  uni¬ 
form  quantiser.  Here  we  make  the  assumption  that  our  quantizer  structure 
is  symmetric  with  respect  to  the  origin.  But  the  optimum  uniform  quantizer 
still  has  its  own  advantages.  It  is  easily  calculable  for  the  x^  (deci¬ 
sion  levels)  and  y^  (quantization  levels).  Also  Wood  [72]  showed  that  it 
gives  the  minimum  distortion  among  the  three  quantizer  structures  if  a 
fixed  entropy  of  the  quantizer  output  is  required. 

Since  we  have  presented  many  tables  for  optimum  nonuniforn,  optimum 
uniform  and  approximate  nonuniform  quantizers,  all  these  tables  are 
created  under  the  assumption  that  we  have  a  unity  variance  input  signal. 
Using  f(n)  to  denote  M.S.E.,  there  is  no  closed-form  expression  to  relate 
f(n)  and  n  (nrnumber  of  quantizer  bits).  However,  one  still  can  easily 
approximate  the  numerical  results  of  f(n)  by  a  piecewise  exponential 
fitting  numerical  technique.  Table  A. 91  shows  the  analytic  models  of 
f(n).  If  we  denote  h(x)  -  f'  *(x),  the  corresponding  analytic  models 
of  h(x)  are  shown  in  Table  A.9II.  Tables  A. 91  and  A. 911  serve  as  a 
useful  tool  for  solving  the  Rate-Distortion  problem  in  Information  theory. 


A. 4  Conclusions 


We  have  obtained  the  approximate  nonuniform  quantizer  characteristics 
for  Causslan  and  I.aplacian  distributed  signals  and  Investigated  the  per¬ 
formance  of  this  quantizer  with  optimum  nonuniform  quantizer  as  well  as 
optimum  uniform  quantizer.  We  have  also  extended  the  performances  of 
optimum  nonuniform  and  optimum  uniform  quantizers  to  large  quantisation 
levels  N  ■  512  (N  ■  2°:n  «*  9)  accompanied  by  useful  information,  eiz., 
mean  square  error,  signal  to  noise  ratio  and  entropy,  etc. 

In  practical  implementation,  approximate  nonuniform  quantizer  seems 
as  good  as  an  optimum  nonuniform  quantizer,  and  it  is  easily  calculable 
with  only  a  marginal  increase  in  its  complexity  compared  to  the  optimum 
uniform  quantizer.  Also  it  is  an  acceptable  and  good  quantizer  structure 
if  the  computer  system  doesn't  have  enough  buffer  space  to  store  the 
entire  tabulated  data  for  every  large  N  or  m. 

The  scheme  described  in  this  Appendix  can  be  easily  extended  to 
other  input  signal  processes  and  is  not  restricted  to  just  Causslan  and 
Laplacian  processes. 
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figure  A  J  dr.n  Square  Error  vereue  Number  of  Quantltar  Bits  for  Signal  with  Cauealan 
Density  (Mean  •  0.  Varlanca  ■  l) 
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figure  A  6  Mean  Square  Error  vereue  Number  of  Quantlier  Blta  for  Signal  with  Laplaclan 


Denelty  (Mean  -  0t  Varlanca  *  l) 
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Figure  A. 7  Signal  to  Noise  Ratio  versus  Number  of  Quantiser  Hits  for  Signal  with  Gaussian 
Penalty  (Mean  -  0,  Variance  -  1) 
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Figure  A  8  Signal  to  Hois#  Rot  io  versus  Niasher  of  Quantizer  Mrs  for  Signal  with  Laplacian 
Density  (Mean  -  0.  V*»r lance  •  1) 
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Figure  A. 11  Quantizer  Output  Level*  and  their  Probability  Dlatrlbutlon  (S 


Table  A.l  Approximate  Average  Amount  of  CDC  Cyber  173  Computer 
Time  Required  to  Determine  the  Quantizer  Structure 
with  Gaussian  Distributed  Innut  Signal 


No.  of  Quantizer  Bit 


CDC  Cyber  173  Computer  Time  (Seconds) 


Approximate 

Nonunlform 

(Compandor) 


Optimum 

Uniform 

Opt lmum 
Nonunlform 

10.193 

382.331 

The  times  shown  are  the  cumulative  computation  times  that  are 
required  to  calculate  the  ouantizers  up  to  the  corresponding  value 
of  n. 


Table  A.  2  Approximate  Average  Amount  of  CDC  Cyber  173  Computer 
Time  Required  to  Determine  the  Ouantizer  Structure 
with  Laplaclan  Distributed  Input  Signal 


No.  of  Quantizer  Bit 


CDC  Cyber  173  Computer  Time  (Seconds) 


Optimum 

Uniform 


Opt imura 
Nonnif orm 


Approximate 

Nonunlform 

(Compandor) 


The  times  shown  are  the  cumulative  comoutatlon  times  that  are 
required  to  calculate  the  quantizers  up  to  the  corresponding  value 
of  n. 
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Optimum  Uniform  Quantisers  for  Signal  with  Gaussian  Density 
(Mean  -  0,  Variance  ■  1) 
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Optimum  Nonuniform  Quantizers  for  Signal  with  Gaussian  Density 
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Table  A. 6  Optimum  Nonuniform  Quantizers  for  Signal  with  Laplaclan 
Density  (Mean  *  0,  Variance  -  1) 
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Table  A. 91  Analytic  Models  for  Some  Often  Used  Quantizers  for 
Unity  Variance  Probability  Density  Functions 


P robability 


Mean  Square  Rrror  Model,  f(n) 
n  ■  number  of  bits 


Quant izer 

Dt>nMty  |  2n  -  number  of  levels 

Shannon 

1  Arbitrary, 

Unitary  Variance 

f(n)  -  2-2n  ,  n  >  0 

Opt imum 
Nonuniform 

Gaussian 

N  (0,1) 

7==-  exp(-x2/2) 

/2tT 

f(n)  -  ( 

2~1 . 504 7n  ,  0  <  2°  <  5 

1.5253*2~1,8274n,  5  <  2"  <  36 
2.2573-2_1'9626n,  36  <  2°  <  512 

Opt imum 

Nonuni  form 

I.aplacian 

Unity  Variance 
^  exp(-/2!x, ) 

f(n)  »  i 

(  2-l.l711n  q  ^  ^  ^ 

2.0851 •2"1’76ASn,  5  <  2n  <  36 
[  3.6308-2-1'qr'72n,  36  <  2n  <  512 

Opt imum 
Nonuniform 

Gaussian 

N  (0,1) 

"ir  exp(-x2/2) 

/2tt 

j 

f(n)  *  ( 

’  2-1.5012n  ,  0  <  2n  <  5 

1.2A77.2_1*6883n,  5  <  2n  <  36 
1.5414-2"1'7562",  36  <  2n  <  512 

Opt imum 
Uniform 

Lapl ac i an 

Unity  Variance 

^  exp(-/2|x: ) 

f (n)  -  < 

2~1 . 1619n  ,  0  <  2n<  5 

1.4156-2~1,4518n,  5  <  2n  <  36 
2.1969-2-1,,i<,44n,  36  <  2n  <  512 

Approximate 
Opt imum 

Non uniform 
(Compandor) 

Gaussian 

N  (0.1) 

■7==  cxp(-x2/2) 

/ffr 

f(n)  -  j 

2~1 -4864n  ,  0  <  2n  <  5 

1.5597-2'1  *8239n,  5  <  2n  <  36 
i  0.2677-2'1*357  \  36  <  2n  <  512 

Approximate 
Opt imum 
Nonuniform 
(Compandor) 

Lap  lac ian 

Unity  Variance 

~  cxp(-/2lx|  ) 

f (n)  -  < 

f  2-1.1491n  ,  0  <  2n  <  5 

2.2164-2-1,7710",  5  5  2n  <  36 
[  4.2461-2“1,9928n.  36  <  2°  <  512 
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Table  A.9II  Formulas  for  h(x)  for  Different  Quantizers 


Quantizer 


Shannon 


Opt  Lmuo 
Nonunlform 
(Gaussian 
Density) 


Opt iraum 
Nonuni  fora 
(Laplacian 
Density) 


Opt iraum 
Uniform 
(Gaussian 
Density) 


h(x)  -  f’(„) 


-1 


h(x)  -  -V  loy>2(-x/2  logp2) 


nil  x  <  0 


0.0403  -  0.6646  log,(-x),  -1.0430  <  x  <  -0.1020 


h(x)  -  ■(  0 .  ‘j  1 99  -  0.5472  1or2(-x),  -0.1020  <  x  <  -0.0027 


0.8247  -  0.5195  log2(-x),  -0.0027  <  x  <  -1.48  x  10 


-5 


h(x)  - 


^-0.2569  -  0.8539  log^-x),  -0.8117  <  x  <  -0.1490 
0.7654  -  0.5667  log2(-x),  -0.1490  <  x  <  -0.0044 


[  1.1753  -  0.5109  1or2(-x),  -0.0044  <  x  <  -2.45  x  10 


-5 


h(x) 


Op  t i nun 
Uniform 
(I.aplacian 
Density) 


0.0382  -  0.6661  log2(-x),  -1.0406  <  x  <  -0.0965 
0.3234  -  0.5923  log2(-x),  -0.0965  <  x  <  -0.0035 


(  0.5170  -  0.5694  lop^(-x),  -0.0035  <x<  -3.28  x  10 


-5 


Approximate 
Optimum 
Nonuni  form 
(Gaussian 
Density) 


-0.2688  -  0.8607  1or2(-x),  -0.8054  <  x  <  -0.1377 
h(x)  -  j  0.3516  -  0.6888  1or2(-x),  -0.1377  <  x  <  -0.0080 


{  0.8026  -  0.6272  1or2(-x),  -0.0080  <  x  <  -1.16  x 


10 


-4 


h(x)  -  < 


0.0290  -  0.6728  1or2(-x),  -1.0303  <  x  <  -0.1047 
0.5371  -  0.5483  1or2(-x),  -0.1047  <  x  <  -0.0019 


(-1.4658  -  0.7368  log2(-x),  -0.0019  <  x  <  -5.29  x  10 


-5 


Approximate 
Opt Inun 
Nonuni  form 
(I.aplaeian 
Density) 


!* -0.2857  -  0.8702  log2(-x),  -0.7965  <  x  <  -0.1573 
i (x)  -  -j  0.8154  -  0.5647  log2(-x),  -0.  1573  <  x  <  -0.0046 
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appfndix  b 


Derivation  of  Approximate  Rantlonal  Function  Operator 


The  SDF  of  an  NCI  model  image  is  given  by 


1  -  rot(z  +  z  1  +  z  +  r  *) 
S(zx.z2) - L~4 - - - if" 2- 

[1  -  a(Zj  +  Zj  +  z2  +  z2  ' i 


A 

Let  R(z^,z2>  ■*  SCz^.z^j)  ,  thus 


RfZj.Zj,)  - 


1  -  a(Zj  +  Zj1  +  z2  +  z"1) 

[1  -  rafzj  +  Zj1  +  z2  +  ]** 


(B-l) 


(B-2) 


We  know  from  the  following  axiom  that 


(l-x)"n  -  1  +  nx  +  SfeS 1  x2  +  niO+lpt2!  x3  +  .  .  . 


(B-3) 


provided  x  <  1  exists. 


Then  the  alternative  expression  of  R(z^,z2)  is,  for  simplicity,  let 

A  ,  -1  ,  -1. 

x  -  a(Zj  +  zJ  +  z2  +  z2  ), 


R(* j »^2) 


(1-rx) 


.  ,  2  2 

(1-x)  1 +»jrx +^-y  ~y~+ . 

,  .,  .  1*  3  r2x2  .  A1 

1  +  *srx  +  — 2 - jy-  +  •  *  •  +  “ 

,  2  1.3  r2x3 

-x  -  «,rx  - 

_  !•  3«  5 •_«.•_( 2n-l)_  rnxn^_ 
-n  nf 


n  n 

-  + 
nl 


l-3-5---(2n-3)  rn_1xn 
2n-l  '  (n- 1  )T 
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Substitution  of  n  -  1,2,...,  into  (B-4) ,  and  after  some  algebraic  aanipu 
lations,  different  order  of  R^n^(Zj,z2)  is  easily  obtained  as 

R(l)(zrz2)  -  1  -  (1  -  |)  a  (Zj  +  zj1  +  z2  +  z"1)  (B-5) 

.  -  1 

0  »  -(1  -^)a 

Figure  B.l  Construction  Diagram  of  R^(Zj,z2> 
r(2>(z1,Z2)  "  1  “  li(1  "  r  +  4>r4a2  ~  <l“f>  Q  (Zj  +z‘1  +  z2  +  Z"1) 

-S(l  -  r  +^)ra2(z^+*“2+z2+z”2  +  2 (z^+z^^+Zj  1z2+z11z21)  ] 

(B-5) 

.  -  l  -  4(1  -  r  +  f)r4a2 

2x  0  2x 

x  0  •  0  x  0  -  -(l-|)a 

2x  0  2x  ~ 

x  -  -4j(l-r  +  J)ra* 
x 

(2) 

Figure  B.2  Construction  Diagram  of  R  (z  ,z_) 
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Figure  R.4  Construction  Diagram  of  R'  (z^ ,Zj) 
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Figure  8.5  Construction  Diagram  of 

It  Is  noted  that  when  the  order  of  R^fipXj)  Roes  up,  wore  sample  points 
In  the  spatial  domain  would  he  Involved  In  calculating  the  approximate 
rational  function  operator  Z- 
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APPENDIX  C 


Unique  Inverse  Solution  of  Function  gj(x) 


The  function  g^(x)  is  given  by 


k;u)  - 


(-21oge2)2"’X 

(1  -  p*2‘2V 


(C-l) 


I.et  h(*)  be  the  inverse  of  g’(*)»  For  convenience  we  denote  g^(x)  and 
— 2x 

2  by  z  and  y  respectively,  and  rewrite  (C-l)  as 


(-21og^2) *y 


i  *  2  .  4  2 

1  -  2P^y  +  p^ 


(C-2) 


■2x 


Since  y  “  2  ,  we  have  1  *  V  > .  0  for  any  given  x  >  0.  Rearranging  (C-2) 

and  getting  an  explicit  expression  as  a  function  of  y,  we  have 


(P^z)y2  -  2(p2z  -log^2)y  +  z  »  0  . 


(C-3) 


The  roots  of  (C-3)  are 


(-P2z  +  log  2)  ♦  J-2p\r.  log  2  +  (log  2)7 

1  »  1  «*  V 


p\z 


(C-4n) 


and 


■(“P.z  +  log  2)  -  /-2pTz  log  2  +  (log  2)‘ 
_ l  0 _ i  »'  e 


(C-4b) 


The  admissible  root  should  lie  in  the  Interval  (0,1 J.  Since  1  >  p2  >  0, 
and  from  (C-l),  we  can  observe  that 


21ogc2 

0>  i> - =— -  for  V  n  >  0 

(1  -  Pj)2 


(C-5) 


i.e.,  z  is  nonpositive.  Therefore 
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(C-6a) 


/^2pTz  log  2  +  (log  2)’  >  0 
It*  c 

-Pj7.  +  loge2  >  0  (C-6b) 

and  it  is  easy  to  prove  that 

-pjz  +  loge2  >  /^[z  loge2  +  (loKe2)r  .  (C-7) 

With  the  above  conditions  (C-6)  and  (C-7) ,  we  conclude 

y2  ~  y\  t  0  •  (c-8) 

In  a  similar  manner,  using  (C-6)  and  (C-7),  it  is  easy  to  check  that 

i)  yj  >  1  ,  hence  y^  is  inadmissible, 
ii)  1  ^  y^  >  0  ,  y ^  is  the  unique  solution. 
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